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At  the  present  state-of-technology  in  the  computational 
aerodynamics,  the  number  of  grid  points  which  can  be  used  in 
solving  a complex  flow  problem  is  often  limited  by  the 
capability  of  existing  computer  system.  Hence,  the 
generation  of  an  optimum  grid  network  is  crucial  to  the 
accurate  solution  of  a complex  flow  problem. 

A general  adaptive  grid  generation  code  based  on  a 
constrained  variational  principle  is  developed  by  minimizing 
a general  variational  functional  which  consists  of 
functionals  for  measuring  the  smoothness,  orthogonality,  and 
concentration  of  grid.  The  coefficients  (penalty  parameter) 
of  the  functionals  of  orthogonality  and  concentration  can  be 
determined  by  dimensional  analysis;  moreover,  they  can  be 
treated  as  variables  to  enhance  the  corresponding  property 
locally.  The  one-dimensional  variational  principle  has  been 
readily  applied  to  adaptive  boundary  grid  generation  as  well 
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as  two-dimensional  case  in  the  domain,  and  excellent 
consistancy  is  achieved  with  the  same  corresponding  penalty 
parameters.  The  developed  one-dimensional  grid  generator 
can  be  applied  to  any  complex  configuration  of  boundary  such 
as  curvilinear  surface  combined  with  polygonal  shape. 

The  adaptive  grid  generation  technique  developed  with 
the  artificial  clustering  is  successfully  applied  to  the 
computation  of  projectile  inviscid  transonic  flow  problem. 

An  excellent  agreement  has  been  obtained  between  the 
pressure  computed  on  the  adaptive  grid  network  and 
experimental  surface  pressure.  Also,  this  agreement  is 
better  than  the  agreement  obtained  from  the  solution 
computed  on  the  fixed  grid  network  which  is  generated  by 
"GRIDGEN." 
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CHAPTER  I 
INTRODUCTION 


In  recent  years,  the  study  of  boundary-fitted  grid 
systems  has  become  important  and  active  because  of  the 
versatile  applications  of  these  systems  in  engineering  pro- 
blems. These  applications  have  made  the  finite-difference 
method  a very  effective  tool  for  solving  complex  fluid  dyna- 
mics problems.  In  general,  a grid  system  in  the  physical 
space  is  numerically  generated  either  by  solving  a set  of 
partial  differential  equations  with  specified  conditions  or 
by  other  methods  such  as  algebraic  interpolation  from  pre- 
determined boundary  points  and/or  along  chosen  grid  lines. 
For  the  boundary-fitted  grid  system,  a general  curvilinear 
coordinate  is  developed,  such  that  all  boundaries  of  the 
arbitrary  two-dimensional  region  of  interest  coincide  with 
coordinate  curves.  Thus,  finite  difference  expressions  at 
the  boundary  may  be  applied  using  exact  grid  points  on  the 
intersections  of  curvilinear  coordinate  lines  without  the 
need  for  any  interpolation  between  points  of  the  grid.  The 
avoiding  of  interpolation  is  particularly  important  for 
boundaries  with  strong  curvature  or  slope  discontinuities, 
both  of  which  are  common  in  physical  applications.  Like- 
wise, interpolation  between  grid  points  not  coincident  with 
the  boundaries  is  especially  inaccurate  with  a governing 
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system  of  equations  that  produces  large  gradients  in  the 
vicinity  of  the  boundaries,  and  the  character  of  solution 
may  be  significantly  altered  in  such  cases  [1],  The  genera- 
tion of  a boundary-fitted  curvilinear  coordinate  system 
which  needs  to  employ  the  technique  of  coordinate  transform- 
ation between  physical  space  and  computational  space  is, 
therefore,  an  important  part  of  a numerical  algorithm  for 
highly  accurate  solution. 

A good  grid  system  generated  for  fluid  flow  computa- 
tions can  be  justified  from  the  smoothness  of  grid  distri- 
bution, the  minimal  skewness  of  grids  and  the  grid 
resolution  being  adaptive  to  solution  in  the  physical 
space.  It  has  been  shown  that  rapid  changes  of  grid  size 
and  highly  skewed  grids  can  result  in  undesirable  errors 
[2].  Also,  it  is  well  known  in  the  approximate  solution 
that  the  choice  of  high  concentration  of  grid  in  regions 
where  the  solution  gradient  is  very  large  is  crucial  to  the 
accuracy  of  the  numerical  result.  In  fact,  an  improper  grid 
resolution  in  high  gradient  regions  can  be  detrimental  to 
the  solution  accuracy  as  well  as  to  the  convergence  process 
of  a computational  algorithm  for  fluid  flow  problems  [3]. 
Therefore,  the  success  of  a finite  difference  method  in  the 
computational  space  for  solving  complex  fluid  dynamics 
problems  can  depend  strongly  upon  the  generation  of  a good 
adaptive  grid  system  in  the  physical  space. 

In  the  last  few  years,  many  different  adaptive  grid 
generation  methods  have  been  proposed  and  reported  in  a 
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symposium  [4],  a workshop  [5]»  and  a conference  [6].  It 
seems  that  each  proposed  method  has  its  own  advantages  for 
the  special  problem  treated,  but  it  may  not  be  suitable  for 
other  complex  problems.  However,  the  method  based  on  a 
constrained  variational  principle,  proposed  by  Brackbill 
[7] i is  more  general  and  shows  great  promise  for  complex 
problems  in  which  physical  phenomena  such  as  shock  waves, 
flames,  and  viscous  layers  may  cause  extremely  large 
gradient  regions  of  unknown  orientation.  Jin  this  method, 
the  governing  equations  of  the  adaptive  grid  generation  are 
derived  by  minimizing  a general  variational  functional. 

This  functional  consists  of  a functional  for  measuring  the 
smoothness  of  grid  and  two  constrained  functionals,  one  for 
the  orthogonality  of  grids  and  the  other  for  the  concentra- 
tion of  grids.  The  relative  importance  of  the  two  penalty 
parameters  on  adaptive  grid  generation  had  been  investigated 
for  simple  problems  and  reported  in  [7],  and  these  two 
parameters  were  treated  as  constants  to  date.  Saltzman  and 
Brackbill  [8,9]  have  applied  the  method  to  a two-dimensional 
inviscid  supersonic  flow  past  a step  in  a wind  tunnel  with 
rectilinear  boundaries;  they  obtained  striking  results  with 
an  adaptive  grid  generator  that  caused  the  computational 
grid  to  move  with  the  shock  fronts.  The  grid  points  on  the 
rectilinear  boundary  are  obtained  by  extrapolating  from  the 
corresponding  interior  grid  lines  in  a simple  linear  rela- 
tion so  that  the  interior  grid  lines  will  be  orthogonal  to 
the  boundary.  This  method  may  be  cumbersome  for  the 
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curvilinear  boundaries  or  more  complicated  body  shapes  such 
as  projectiles. 

A better  design  of  flight  vehicles,  aerodynamic 
devices,  and  sheltering  structures  can  depend  strongly  upon 
an  accurate  prediction  of  aerodynamic  forces.  There  are 
quite  a few  effective  f inite-difference  schemes  and  solution 
algorithms  developed  for  solving  complete  and/or  modified 
Navier-Stokes  equations;  the  choice  of  a good  grid  system 
for  the  approximation  may,  however,  be  crucial  to  the  suc- 
cess of  these  routines.  For  a complex  flow  problem,  an 
extremely  large  number  of  grids  is  often  required  to  attain 
satisfactory  results;  consequently,  an  optimum  or  proper 
adaptive  gridding  becomes  essential  for  accurate  aerodynamic 
force  computations  wqth  limited  computer  resources. 

Steger  [10],  at  NASA  Ames  Research  Center,  developed  a 
thin-layer  Navier-Stokes  code  about  arbitrary  two-dimen- 
sional geometries.  Later,  Pulliam  and  Steger  [11]  extended 
the  thin-layer  code  to  the  three-dimensional  case.  A combi- 
nation of  general  curvilinear  coordinate  transformations  and 
an  implicit  approximate  factorization  technique  which  was 
developed  by  Warming  and  Beam  [12]  was  employed  in  that 
code;  so  that  the  fine  grid  sizes  required  for  spatial 
accuracy  and  viscous  resolution  do  not  impose  stringent 
stability  limitations.  The  present  state  of  the  technology 
in  computational  aerodynamics  clearly  indicates  that  the 
thin-layer  code  can  provide  accurate  solutions  for  high 
speed  aerodynamics  problem  [11,13].  Also,  Nietubicz, 
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Pulliam,  and  Steger  [14]  have  successfully  applied  the  thin- 
layer  Navier-Stokes  code  to  solve  ax i symmetric  flow  pro- 
blems. They  modified  the  code  by  use  of  a cylindrical 
coordinate  system,  and  the  three-dimensional  governing  equa- 
tions are  simplified  for  flow  fields  which  are  invariant  in 
the  circumferential  direction.  Moreover,  this  axisymmetric 
thin-layer  code  has  readily  been  applied  to  transonic  pro- 
jectile aerodynamics  problems  at  the  U.S.  Army  Ballistic 
Research  Laboratory  [15].  Their  report  shows  that  good 
agreement  has  been  obtained  between  the  computed  and  experi- 
mented surface  pressure.  However,  the  successful  applica- 

% 

tion  of  the  code  depends  strongly  upon  a good  adaptive  grid 
network  generated  from  a grid  generation  code  named  GRIDGEN 
[16].  For  the  procedure  of  grid  generation  in  GRIDGEN,  th§ 
boundary  grid  points  are  generated  by  the  clustering  which 
is  implemented  with  a cubic  function  that  allows  the  user  to 
specify  the  first  and  the  last  grid  point  spacing  along  the 
grid  line  of  each  segment  of  the  boundary.  A planar  grid 
network  within  the  domain  can  be  generated  by  solving  either 
an  elliptic  boundary  value  problem  or  a hyperbolic  initial 
value  problem  in  the  current  version  of  GRIDGEN.  The 
resulting  grid  network  can  then  be  modified  to  provide  high 
grid  resolution  in  the  viscous  layer  by  clustering  grid 
points  along  the  grid  lines  normal  to  the  streamwise  direc- 
tion. Finally,  a three-dimensional  grid  network  for  a pro- 
jectile aerodynamic  problem  is  comprised  of  a sequence  of 
the  generated  planar  grid  networks  around  the  axis  of  the 
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projectile.  As  reported  by  Sturck  [17],  an  accurate  solu- 
tion for  the  aerodynamics  of  an  axisymmetric  transonic  pro- 
jectile can  be  obtained  with  a good  adaptive  grid  network. 
This  grid  network  is  generated  by  GRIDGEN  only  after 
considerable  experimentation  with  boundary  grid  positioning 
and  grid  density  distribution. 

For  a transonic  projectile  at  an  angle  of  attack,  the 
three-dimensional  grid  network  generated  by  GRIDGEN  may  not 
be  a proper  one  for  use  in  the  computation,  since  the  flow 
field  is  not  axisymmetric  and  consequently  a planar  grid 
network  adaptive  to  the  flow  field  at  an  azimuthal  location 
may  not  be  a good  grid  network  for  use  at  other  azimuthal 
planes.  Moreover,  for  projectiles  at  supersonic  speed  and 
other  high  speed  viscous  flow  problems  such  as  blast  aero- 
dynamics [18],  the  regions  where  the  solution  gradients  are 
very  large  are  not  known  in  advance;  consequently,  a cut  and 
try  approach  for  finding  a good  grid  network  can  be  very 
tedious.  Therefore,  it  is  extremely  important  and  essential 
for  solving  complex  flow  problems  to  develop  effective  tech- 
niques for  systematic  generation  of  good  adaptive  grid  net- 
works. 

The  prime  purpose  of  this  study  is  to  investigate  and 
develop  a two-dimensional  adaptive  grid  generation  technique 
for  high  speed  projectile  aerodynamics  computation.  The 
secant-ogive-cylinder-boattail  projectile  with  sting  and  the 
axisymmetric  version  of  the  thin-layer  Navier-Stokes  code  of 
the  U.S.  Army  Ballistic  Research  Laboratory  are  employed  in 
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this  development.  The  adaptive  grid  system  is  mainly  based 
on  a constrained  variational  principle;  hence,  the  governing 
equations  for  the  transformation  are  quasi-linear  and  must 
be  solved  numerically  by  a stable  iterative  method.  The 
detailed  investigation  is  performed  in  this  study  for 
several  different  solution  methods.  The  method  of  Newton- 
Raphson  iteration  with  successive  over  relaxation,  which  has 
been  shown  to  be  a highly  efficient  and  accurate  method  in 
solving  nonlinear  boundary  value  problems,  is  employed  here 
to  solve  the  set  of  complicated  partial  differential  equa- 
tions. The  investigations  conducted  on  the  constant  penalty 
parameters  for  orthogonality  and  concentration  of  grid 
points  have  shown  that  these  two  parameters  can  be  treated 
as  variables  over  the  field  by  using  dimensional  analysis  to 
enhance  the  properties  of  the  control  function  obtained  from 
the  approximate  solutions.  The  adaptive  boundary  grid  gene- 
ration developed  by  applying  a one-dimensional  variational 
principle  provides  complete  consistency  between  the  grid 
points  within  the  domain  and  corresponding  points  on  the 
boundary.  The  developed  adaptive  grid  generation  code 
coupled  with  an  axisymmetric  thin-layer  Navier-Stokes  code 
has  demonstrated  the  need  for  clustering  the  grid  points  in 
the  vicinity  of  the  body  surface  to  obtain  good  results. 

This  coupled  code  gives  computational  results  which  agree 
well  with  the  experimental  data. 


CHAPTER  II 

ADAPTIVE  GRID  GENERATION 


II . 1 Two-Dimensional  Coordinate  Transformation 
For  a two-dimensional  coordinate  system,  the  general 
transformation  can  be  an  one-to-one  mapping  of  an  irregular 
boundary  curve  surface  in  the  physical  plane  (x,z)  onto  a 
rectangular  boundary  in  the  computational  plane  (?,?),  as 
shown  in  Fig.  2.1.  The  solution  in  the  transformed 
computational  plane  is  then  given  by  the  functions 


5 = 5(x,z) 
C = c(x,z) 


(2.1) 


The  grid  of  the  computational  plane  is  a square  mesh  with 
equal  grid  spacing  in  a fixed  rectangular  field  with  MxN 
grid  points.  The  fluid  flow  problems  are  then  solved  on 
this  computational  plane.  Similarly,  the  inverse 
transformation  is  given  by 


x = x( 5, e) 

z = z(c,c) 


(2.2) 


where  5,  5 are  the  natural  coordinates.  The  interior  grid 
points  in  the  physical  plane  are  the  intersection  points  of 
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Fig.  2.1  Coordinate  transformation 
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the  natural  coordinate  curves  g = const,  and  5 = const., 
which  intersect  the  boundary  in  the  chosen  boundary  grid 
points.  This  forms  the  boundary-fitted  grid  system.  Thus, 
the  integer  values  j and  k of  the  natural  coordinates  (5,5) 
define  a point  (j,k)  in  the  computational  plane,  and  the 
mapping  associates  with  it  a point  in  the 

physical  plane.  In  formulating  the  grid  generation  problem, 
mathematically  the  grid  [x^k,Zj^]  is  viewed  as  the  image 
of  a mapping  [ x( £ , ? ) , z( £ , ? ) ] in  which  only  the  points 
corresponding  to  the  integer  values  of  the  natural 
coordinates,  (£,£)>  are  realized.  Conversely,  the  image  of 
a computation  grid  of  curvilinear  quadrilateral  cells  in  the 
physical  plane  is  a uniform,  rectilinear  grid  in  the  (5,5) 
plane  with  spacing  defined  as 

AC  = Ac  = 1 (2.3) 

A grid  generator  determines  the  mapping  [x( £ , 5 ) , z( £ , 5 ) ] . 

The  system  of  governing  equations  for  the  mapping  can 
be  defined  by  numerical  transformation  procedures  including 
techniques  based  on  elliptic  and  hyperbolic  partial 
differential  equations,  algebraic  or  geometric  procedures 
based  on  parametric  surfaces  or  interpolation  techniques, 
and  analytic  transformation  procedures  such  as  conformal 
mapping.  To  generate  an  adaptive  grid  system,  in  this  work 
we  use  partial  differential  equations  as  the  governing 
equations  for  the  transformation  of  the  physical  plane  to 
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the  computational  plane.  We  then  define 


L[ g(x,z) ] = 0 
L[c(x,z)]  = 0 


(2.4) 


where  L[  ] is  an  operator;  thus  Eqn.  (2.4)  represents  some 
type  of  partial  differential  equations  (PDE)  with  Dirichlet 
boundary  conditions  such  that  ? = 1 , ? = N,  g = 1,  and  g = M 
of  the  computational  plane  are  on  the  corresponding  segments 
AE,  BD,  AB,  and  DE  of  the  physical  plane  (see  Fig.  2.1). 
Since  it  is  desired  to  perform  all  numerical  computations  in 
the  uniform  rectangular  transformed  plane,  the  dependent  and 
independent  variables  in  Eqn.  (2.4)  must  be  interchanged . 
This  yields  the  coupled  system  of  equations  in  the  L*[  ] 
operator  form 


L*[x(s,?)]  = 0 
L*[z(g,5)]  = 0 


(2.5) 


These  equations  are  subject  to  the  transformed  boundary 
conditions  which  are  suitably  chosen  functions  that  map  the 
known  shape  of  four  segments  AE,  AB,  BE,  and  DE  of  the 
physical  boundaries  into  the  corresponding  segments  AE  ( 5 = 
1),  AB  (5  = 1),  BD  (5  = N),  and  DE  (5  = M)  on  the 
rectilinear  boundaries  of  the  chosen  computational  domain. 

Usually,  the  system  given  by  the  equations  of  (2.5)  is 
considerably  more  complex  than  that  defined  by  Eqn.  (2.4); 
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however,  the  boundary  conditions  of  Eqn.  (2.5)  are  specified 
on  straight  boundaries  and  the  coordinate  spacing  in  the 
computational  plane  is  uniform. 

II .2  Review  of  Variational  Principles 
A variational  principle  specifies  a scalar  functional, 
I*,  which  is  defined  by  an  integral  form 

I*  = /fiA(U)da  + /rB(U)dr  (2.6) 

where  A and  B are  the  operators  such  that  A(U)  and  B(U) 
represent  the  constraint  equations  in  terms  of  the  unknown 
vector  function  U (i.e.,  U = {U}),  in  a domain  n,  and  on  the 
boundary  r.  The  solution  to  the  continuum  problem  is  a 
function  U which  makes  I*  stationary  with  respect  to  small 
changes  6U,  or  in  other  words,  the  function  U which 
minimizes  I*  is  in  some  sense  the  "best"  choice.  Thus,  for 
a solution  to  the  continuum  problem,  the  variation  is 

6 1*  = 0 

If  a variational  principle  can  be  found,  then  immediately 
means  are  established  for  obtaining  approximate  solutions  in 
the  standard,  integral  form  suitable  for  adaptive  grid 
generation  analysis. 

A variational  principle  may  also  be  considered  as  the 
problem  of  making  a functional  I*  stationary,  subject  to  the 
unknown  vector  U obeying  some  set  of  additional  differential 
constraint  equations 
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C(U)  =0  in  ft 
E(U)  = 0 on  r 


(2.7) 


where  C and  E are  the  operators  such  that  C(U)  = { C ( U) } and 
E(U)  = |E(U)}.  We  can  introduce  these  constraint  equations 
by  forming  another  functional 


1 = X*  + U cT(u)c(u)d«  + *r  Jr  ET(U)E(U)dr  (2.8) 

in  which  x^  and  xr  are  the  penalty  parameters  or  so-called 
"penalty  numbers"  in  reference  [19],  and  superscript  "T" 
stands  for  transpose  of  the  matrix.  The  products,  C^C  and 
EtE,  are 


CTC  = [c^,c|,  . . .] 

ETE  = [E^,E^,  . . .] 

and  the  quantities  of  CT(U)C(U)  and  ET(U)E(U)  must  always  be 
positive.  Substituting  Eqn.  (2.6)  into  Eqn.  (2.8)  gives 

1 - In  * *r  (2-9) 

where  In  is  the  functional  within  the  domain  n, 


in  = /nA(U)dn  + /ncT(U)C(U)dn, 


(2.10) 
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and  Ip  is  the  functional  on  the  boundary  r, 


Ir  = / pB(U)dr  + x r /rET(U)E(U)dr, 


(2.11) 


II .3  Governing  Equations 
Since  f and  ? are  curvilinear  coordinates  in  the 
physical  plane  as  shown  in  Fig.  (2.1),  the  differential 


can  be  defined  from  the  relationship  between  the  Cartesian 
and  curvilinear  coordinates,  and  they  are  directed 
tangentially  to  the  c = const,  and  £ = const,  curves, 
respectively.  A useful  observation  is  that  the  differential 
properties  of  the  mapping  determine  the  properties  of  the 
computations  grid.  For  example,  the  magnitude  of  the  vector 
{D5}  can  be  related  to  grid  spacing  As  by  taking  the  finite 
difference  form  of  the  Eqn.  (2.12)  with  A?  = A?  = 1 , 


vectors 


(2.12) 


[(Ax)2  + (AZ)2]^2  = As 

J 9 & 


(2.13) 


where  As  is  the  grid  spacing  in  the  physical  plane  between 
two  points  (j,k)  and  (j+1,k).  Similarly,  since  the  length 


of  the  vector  resulting  from  a cross  product  of  {Dc}»{D?}  is 
related  to  the  area  of  the  elementary  parallelgram  in  the 
curvilinear  coordinate  system,  we  obtain 


d(Area)  = J~^d£d? 


(2.14A) 


where  J“^  is  the  Jacobian  determinant,  defined  as 


(2.15) 


Thus,  the  finite  difference  form  of  Eqn.  (2.14A)  becomes 


As  described  in  Chapter  I,  a good  grid  system  generated 
for  fluid  flow  computations  can  be  justified  from  the 
smoothness  of  grid  distribution,  the  orthogonality  of  grids 
and  the  grid  resolution  being  adaptive  to  a control  weight 
function  obtained  from  the  solution  in  the  physical  space. 
Thus,  the  governing  equations  of  the  adaptive  grid 
generation  are  derived  from  minimizing  a general  variational 
functional  as  defined  in  Eqn.  (2.10) 


A(Area)  = = j 


(2.14B) 


grid  area  in  physical  space 
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The  functional  for  measuring  the  smoothness  of  grids  is 
Is  which  can  be  written  as 

Is  = fn  C ( ve • ve ) + (vc-vc)]dxdz  (2.17) 

Thus,  minimizing  this  integral  results  in  Laplace's 
equations  for  ? and  ? as  functions  of  x and  z.  Because  of 
the  well-known  averaging  properties  of  solutions  of 
Laplace's  equation,  we  might  expect  a grid  constructed  in 
this  way  to  be,  in  some  sense,  smooth  [20]. 

The  functional  for  measuring  the  orthogonality  of  grids 

is  IQ 


!0  = Jn  ( v?-ve)I 2(y)3dxdz  , (2.18) 

in  which  (V£»v?)  is  zero  when  the  conjugate  curves  (C  = 
const.,  c = const.)  of  the  grid  are  completely  orthogonal. 
Also,  the  inclusion  of  the  (— ) , the  cubic  of  the  Jacobian, 
weights  the  measure  in  favor  of  orthogonalizing  regions  of 
large  cell  area  [9];  this  weighting  function  can  also  reduce 
problems  with  rounding  errors  [7]. 

The  functional  for  measuring  the  concentration  of  grids 

is  Iw 


I = /_  WJ-1dxdz 

W U* 


(2.19) 


where  W = W(x,z) 


is  a given  weighting  function  and 


represents  the  grid  area  in  the  curvilinear  coordinate 
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system  as  derived  in  Eqn.  (2.14B).  If  we  were  to  minimize 
this  integral,  we  would  predict  that  where  W is  large,  J-"' 

A 

(cell  area)  should  be  small  and  conversely  where  J is 
large,  W should  be  relatively  small.  The  parameters,  A0  and 
Aw,  are  specified  coefficients  whose  values  give  more  or 
less  emphasis  to  orthogonality  or  to  concentration, 
respectively.  The  procedure  for  optimizing  the  grid  is  to 
choose  a mapping  that  minimizes  the  functional  In  of  Eqn. 
(2.16) . 

It  is  useful  to  normalize  the  functionals,  IQ  and  Iw, 
in  order  to  keep  the  order  of  magnitude  of  each  functional 
the  same.  Thus,  three  normalization  factors,  \'Q,  A^,  and  A 
are  defined  by 


in  which  \ is  the  specified  global  constant  parameter.  A 
fairly  effective  and  simple  way  to  determine  \'Q  and  A^  is  to 
apply  a dimensional  analysis.  For  instance,  the  derivatives 
. . . , etc.  are  proportional  to  Vl',  where  L is  the 
physical  length  and  L'  is  the  parameter  length  (points 
number).  The  order  of  magnitude  of  the  Jacobian  can  be 
written  as 


9 z 9x  / L x 2 
35  35  ~ 'L,  ; 


9 X 9 Z 

35  35 


(2.21 ) 
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The  integrand  of  the  smoothness  functional  in  Eqn.  (2.17) 
has  the  order  of  (L'/L)^  and,  similarly,  the  integrands  of 
the  orthogonality  and  concentration  functionals  are  (L/L’)^, 

— o — 

¥(L/L')  » respectively,  where  ¥ is  the  unit  of  measure  of 
the  weight  function  within  the  domain.  Thus,  the  order  of 
magnitude  of  the  functionals  is 


Is  ~ (^-)2L2  = (L')2 


I ~ ¥ (— — ) 2L2 
w L* 


(2.22) 


(2.23) 


In  order  that  IQ  and  Iw  have  the  same  order  of  magnitude  as 
Is,  and  A^  should  have 


A ' ~ (ii-)  4 

0 L* 

X.  . ¥ (— — ) 4 

w 

W L* 


(2.24) 


Specific  choices  for  L,  L' , and  ¥ were  made  by  Saltzman 
and  Brackbill  [8];  L was  the  maximum  length  in  the  physical 
domain,  L*  was  the  number  of  points  in  one  direction,  and  ¥ 
was  the  area  average  of  ¥(x,z),  such  that  ¥ = /fi¥ dft. 

Thus,  Eqn.  (2.1b)  becomes 
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I 

ft 


zs  + 


xo 


+ I 


xw 


w 


= fQ  (v?-V5  + VC 


+ — J„wJ'1da  (2.25) 

X ' n 
Aw 

Next,  the  independent  and  dependent  variables  are 
interchanged  by  the  relations 

‘ 1 0 
0 1 


x x 
X z 


z z 
X z 


X X 

5 C 


z^  z 
.5  c J 


?x  ^z 


«X  ?Zj 


(2.26) 


where  , cx  = , . . . , etc.,  are  the  metrics. 

Thus,  they  have  the  following  relations 


^x 

«x 


(2.27) 


where  J is  the  Jacobian  of  transformation  such  that  J = 
in  which  J-1  is  defined  in  Eqn.  (2.15).  The 
functionals  of  Eqn.  (2.25)  can  be  readily  derived  by  using 
the  relations  of  Eqn.  (2.27);  the  variational  integral 
becomes 
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I,  = / , J( vx. vx+vz. vz)dq'  + 

a 6 it 

f . (x  X +2  z )2d n'  + 

Ao 

f , W(x,z)(^)2dq'  (2.28) 

Aw 

where  da'  = dgd^  is  the  area  element  in  the  computational 
domain.  The  necessary  conditions  for  x(5,s)  and  z(5,s)  to 
minimize  I^i  give  the  Euler  equations 


3 3_ 

35  3x? 


Jj— ] [J(  vx-vx+vz- vz)  + 


77  ‘WVP' 

o 


)] 


0 = [- - - 

3Z  35  3z 


3 5 3Z 


-]  [ J(  vx . vx+ vz  • vz)  + 


( x x + z z )2  + ^—  (^4r)] 


‘i 


5 5 5 5 


Performing  the  above  differentiations  and  collecting 
coefficents  of  the  highest  derivative  terms  result  in  the 
coupled  quasi-linear  elliptic  system  of  equations 
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a1X5C+a2xU+a3X^+b1Z5C+b2z^+V^ 


-(— )- 


W 


X'  2J 
w 


b1X5C+b2XU+b3X^+C1Z?5+C2Z?c+C3Z^ 


-(— )- 


W 


X'  2J 
w 


(2.29) 


where  W = + W c and  W = W_?  + can  be 

x £ x ? x z ? z ? z 

transformed  by  applying  Eqn.  (2.27)  to  the  computational 
form 


Wx  - J(Z;W5  - z5W£) 


Wz  = J(x5W£  - x£W?) 


(2.30 


The  coefficients  in  Eqn.  (2.29)  are  given  by 


i = 1,2,3 


(2.31) 


in  which  the  corresponding  coefficients  of  smoothness, 
orthogonality  and  concentration  are 


as1 

as2 

a - 

s3 

aA* 

-2BA* 

YA*  ' 

bs1 

bs2 

b , 
s3 

= 

-cxB* 

2BB* 

-yb* 

.CS1 

cs2 

cs3  _ 

aC* 

-2BC* 

YC*  _ 

* 


(2.32) 
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"aoi 

ao2 

ao3 

2(2x5xc+z5zc) 

bo1 

bo2 

bo3 

= 

Vc 

X5ZC+XCZ5 

V? 

* 

_co1 

co2 

Co3 

z2 

c 

2(xsx;+2z5z5) 

Zf 

wl 

aw2 

aw3 

2 

zc 

-2z5z; 

2 

ZC 

wl 

bw2 

bw3 

= W 

-x  z 

c ? 

XsVVs 

"xcz5 

> 

wl 

Cw2 

cw3 

C\J  VJ> 
X 

1 

-2Vc 

x2 

5 

respectively,  and 

A*  = + z^  , B*  = x^z^  + x?z?  , C*  = x^  + x^ 

a = (x^  + z^)/J5  , 8 = (x^x?  + z5z?)/J3  , 

Y = (x2  + Z2)/j3 


(2.33) 


(2.34) 


(2.35) 


For  convenience,  Eqn.  (2.29)  may  also  be  written  in  the 
operator  form, 


W i ? 

LA[x]  + LB[z]  = -Aw(-|)(j)2 

LB[x]  + LC[z]  = -Aw(-|)(-L)2 

where  X = X/X'  , and 
w w 


(2.36) 
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T . _ 3 a a 

LA  = a„  — 7?  + an  + a 


1 3£2  2 3£3?  ^ 3 £2 


LB  = b„  — — + b 


1 352  2 3C3C  2 9c2 


T„  _ a a . a 

LC  - c.  o + Co  + c.,  q 

1 35^  ^ a c a c 5 a^ 


(2.37) 


II .4  Boundary  Conditions 
The  consistency  between  boundary  and  domain  is 
important  for  the  success  of  the  adaptive  grid  generation. 
Thus,  it  is  necessary  to  generate  the  adaptive  grids  on  the 
boundary.  Saltzman  and  Brackbill  [8]  have  applied  the 

method  based  on  a variational  principle  to  a two-dimensional 

* 

inviscid  supersonic  flow  past  a step  in  a wind  tunnel  with 
rectilinear  boundaries  and  obtained  striking  results. 
Boundary  grids  are  extrapolated  from  interior  grids  of  the 
domain  orthogonal  to  the  boundary.  However,  this  is 
difficult  to  apply  to  extrapolating  on  a curvilinear 
boundary  such  as  a projectile  with  curved  ogive  or  some 
other  irregular  body  configuration.  Therefore,  a one- 
dimensional variational  procedure  must  be  applied  to 
properly  determine  the  adaptive  boundary  grid,  as  follows. 

Similarly  to  the  coordinate  transformation  as 
described  in  sec.  II. 3,  let  S be  a variable  along  a curved 
boundary  in  physical  space.  Then,  S is  mapped 
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correspondingly  to  a variable  5 on  a line  c = const,  in  the 
computational  space  as  shown  in  Fig.  2.2,  and  conversely, 
such  that 


£ = 5(S)  = j 

S = S(5)  = S.  (2.38) 

J 

<3  = 1,2,  • • • , M 

The  first  derivative  Sg  can  be  considered  a measure  of  the 
grid  spacing  AS  by  the  mean  value  theorem,  since  A?  = 1 . 

43  - Sj  - So-1  ■ ^df1  - Sj-1  < s*  < s3  <2-39> 

Thus  represents  the  grid  spacing  in  one-dimensional 
space,  analogous  to  J-"'  , the  cell  area  in  two-dimensional 
space  for  A£  = Ac  = 1 . 

To  optimize  the  physical  characteristics  of  grid 
smoothness  and  concentration  on  the  boundary,  we  seek  to 
minimize  their  functionals  in  accordance  with  the  one- 
dimensional variation  principles.  Thus,  the  functional  of 
Eqn.  (2.11),  which  measures  the  properties  of  the  mapping, 
becomes 


I r = I + x 1 
s w w 


(2.40) 


where  I . the  functional  of  smoothness,  is 
s 


i - /IM  UJ2  ds 

s 1 s^  s 


(2.41) 
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Z 


£ 


Fig.  2.2  1-D  Coordinate  transformation 
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and  I , the  functional  of  concentration,  is 
w 


I 


w 


Lm  V/(S)S  ds 
S1  5 


(2.42) 


in  which  W(S)  is  a specified  weight  function  which  makes  the 
boundary  grid  spacing  small  corresponding  to  the  large  value 
of  W(3).  The  parameter  Aw  also  becomes  A = A/A^  as  defined 
and  determined  in  the  two-dimensional  case  by  the 
dimensional  analysis. 

The  order  of  magnitude  for  the  functionals  of 
smoothness  and  concentration  are 


(2.43) 

(2.44) 

where  L,  a total  length  of  boundary,  can  be  determined  by 
L = fT  dr,  L'  is  the  number  of  corresponding  points,  and  ¥ 
is  the  length  average  of  W(S),  such  that  W = — fr  Wdr. 

Thus,  A^  is  determined  by  matching  the  order  of  Ig  to  1^ 

K = W(rr)3  (2-45) 

W L 

After  interchange  of  the  independent  and  dependent 


(£->2l 


I - W(— )L 


w 


variables,  Eqn.  (2.40)  becomes 
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Ke)  = /?  (S J"1  d5  + — W* ( S ) 2 d?  (2.46) 

'5  \ t £ 

W 

The  necessary  condition  for  S($)  to  minimize  1(5)  gives 
Euler's  equation 


0 = 


[— 

L 3S 


_ 2 2_i  r Q 

35  as.JLb 


-1 


— (WSJ)] 

A ' 4 

w 


(2.47) 


Collecting  coefficients  of  the  highest  derivative  yields  the 
equation, 


S„  = -(— ) W _sf / (2(1  + — S^W)]  (2.48) 

xt  s 5 x,  5 

w w 

where  Ws  = W^/S^,  and  the  end  points  are  5(S^)  = 1,  ?(SM)  = 
M.  Notice  that  A may  be  the  same  as  that  used  in  Eqn. 
(2.29),  because  of  the  requirement  of  consistency  between 
boundary  and  domain. 

II .5  Methods  of  Solution 

In  order  to  generate  a grid  network,  the  governing 
equations,  Eqn.  (2.36),  are  solved  by  the  finite  difference 
method.  The  values  of  the  derivatives  at  the  nodes  are 
approximated  by  the  second-order  central  differences.  For 
the  use  of  an  alternating  direction  implicit  method,  Eqn. 
(2.36)  can  be  approximated  in  either  the  5-  or  5-direction, 
such  that  for  the  k line  (row): 
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a1Xj-1 ,k-2(a1 +a3)x j ,k+a1 Xj+1 ,k+b1zj-1 ,k"2(b1+b2)zj,k 
+b1zj+1,k  = <Sx>d,k 

b1xj-1 ,k'2(b1+b3)xj,k+b1xj+1 ,k+c1zj-1 ,k~2(c1 +c2)zj ,k 
+c1zj+1 ,k  - ^Sz^j,k 


(2.49) 


where 


(Sv^  v = ~ ~ +b0z  1 


x j,k 


i 2 j x j j , k Lr5cr5cjj,k 


w 


" a3[xj,k-1+xj,k+1]  " b3[zj,k-1+zj,k+1] ’ 


= - T,”[b0x  +c0z,._]  . 


z j,k 


. 2VJ/  zJ,j,k  LU2A{;c^2^eJj,k 


w 


b3[xj,k-1+xj,k+1]  " c3[zd,k-1+zo,k+1]» 


J.T- 

and  for  the  j line  (column) 


'3xdA-1-2(a1ta3)xa,k+Vd,k+1tb3zd,k-1-2(b1tb3)zj,k 


+b3zj jk+1  - (si)j,k 


(2.50) 


b3Xa,k-l"2(b1+b3)xj,k+b3Xj,k+1+c3Z,j,k-l‘2(c1+c3)zj,k 
+C3Zj,k+1  = (Sz}j,k 


where 
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<si>0,k 


Aw 


a1  [xj-1 ,k+xj+1 ,k]  " b1 [zj-1 ,k+za+1 ,k] ’ 


(S' ) . . 
z'j,k 


77[i(T)2wz]j,k-tb2x^+c2zU]J,k 

w 


" b1[xo-1,k+xj+1,k]  " 


°1 [zd-1 ,k+zj+1 ,k] 


in  which  ¥x  and  Wz,  the  derivatives  of  the  weight  function 
are  also  approximated  by  the  2nd  order  central  difference, 


‘Vj.k  - [J(zcVZ5Wp]j,k 


Jd.k^zd.k+1”zd,k-1  ^Wj+1  ,k_v,j-1  ,k^ 

(zj+1 ,k"zj-1 ,k)(W3,k+1"Wo,k-1)'/4 


(Vj,k  ' [J(x5Wc-x5W5)].jk 


J;))k[(xJ+1,k-Xd-1,k)(Wj,k+1-Wd>k-1)  - 


(xJ,k+1'xj,k-1  ,(wj+1  ,k'WJ-1,k):l/4 


Also,  the  coefficients  of  Eqn.  (2.31)  which  are  in  terms  of 
the  derivatives  are  approximated  by  the  central  difference 
form. 

There  are  many  methods  to  solve  the  Eqn.  (2.49)  or  Eqn. 
(2.50).  Since  the  governing  equations  are  a coupled 
quasilinear  set  of  equations,  the  finite  difference 
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approximations  are  solved  by  iteration.  In  order  to  obtain 
an  effective  calculation  for  solving  the  grid  system 
equation,  we  investigate  several  methods  to  identify  which 
one  is  the  best.  Methods  of  Newton-Raphson  iteration, 
implicit  line  iteration,  and  alternating  direction  implicit 
are  the  quite  general  and  efficient  solvers.  These  methods 
continually  iterate  with  successive  overrelaxation  (SOR) 
until  the  residual  errors  of  the  solutions  are  everywhere 
less  than  the  allowed  error. 

II. 5.1  Newton-Raphson  Iteration  (N-R) 

Let  Rx  and  Rz  be  the  residual  errors  of  the  Eqn. 

(2.36),  that  is, 

(Rx)d>k  = fLA[x]  ♦ LB[z]  ♦ (i-)|£(l)2}.>k 

A 

w 

(2.51) 

(Rz)^,  = {LB[x]  ♦ LC[z]  ♦ 

w 

In  the  iteration  process,  the  values  of  Xj^  and  zj^  are 
calculated  point  by  point  from  Eqn.  (2.51),  such  that  first 
we  consider  a point  (x^,zA)j  ^ which  is  not  a root  of 
equations  of  (Rx)j>k  = 0 and  (Rz)^^  = 0,  but  is  reasonably 
close  to  the  roots,  where  the  superscript  l is  the  iteration 
number.  We  expand  RxA+1  and  Rzi+1  in  a Taylor  series  about 
xj  k and  Zj,k»  anc*  retain  only  the  linear  terms. 
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(Rx)‘+,1  ■ (R*)i  t+(x‘T,-x*K  „[ 


£ + 1 £■ 


3 (Rx)  -i  k n 

J ? *- 


j,k 


d*k 


j,kL  3X.>k 


(z£+1_z£)  [ 3 & 
( ;j,kL  azdfk  J 


(Rz)j,k  (Rz)j,k+u  -x  }a,kL  3x.jk  J 


(2.52) 


+ (z£+1-z£)  9 ( Rz ^ j | k ] £ 

( }j,kL  3zJ>k 


where  (Rx)*>k  = Rx(x^k,zJ>k)  and  (Rz)^k  = Rz(xJfk,zJfk) , 
in  which  j = 2,  . . M-1  and  k = 2,  . . N-1 . An 

approximate  value  of  the  root  (x!'+\zll'+^  ) ^ k can  be  obtained 
by  setting  Rx£+1=  Rz£+1  = 0 at  the  left-hand  side  of  (2.52) 
to  yield 


o - (rx)£  +(x£+1-x£)  r8(Rx).ii-k-]* 

U - Uxj.)kHx  x ).fki  3x^k  J 


+ (z*+1-z£)  [8(RX--j-xii]  1 

( ^j.k1  3Z.>k  J 


. 3 ( Rz ) . , 

0 = (Rz)  ^ k+(x£+1-x*),  k[ 

J,k  J,K  3Xj}k 

3 (Rz)  . 


+ U*+1- z‘),  J 


J’k  3zo,k 


ilk.j£ 


(2.53) 


where  the  four  derivative  terms  can  be  readily  found  by 
dif ferentiating  the  finite  difference  form  of  Eqn.  (2.49)  or 
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(2.50)  but  note  that  there  are  no  contributions  when 
differentiating  their  source  terms,  (SY).  w (s„^  ,,  or 

ZJjK 

^sx^,k»  ^Sz^j,k  aCCsx) j,k^/axj,k  = 3t(sx^ j,k^/3zj,k 

= . . . = 0) . Thus,  there  results 


9 (Rx)  . , 

g-  = -2(a  +a  ) 

Xj,k  1 3 


k 


a(Rx).i,k  3(Rz)- 
&z  ‘ = 


- -2(b1+b3>J.l 


(2.54) 


3(Rz)i  k 

-Tz  3>  - -2(c  +c  ) . 

zj,k  1 3 


£ ,*1  £ + 1 

Solving  for  x^k  and  z.  ^ from  Eqn.  (2.53)  gives 


xUl 

j,k 


zUl 

0,k 


A 1 3(Rx)i  k 9 (Rz)  . , j. 


(2.55) 


l 

z . , + 


3-k  Dj,k 


1 3 (Rz)  . t 3(Rx)  . , , 

t(Rx),  -(Rz),  v-J,k] 


J’k  8x0,k 


J’k  3xd,k 


where 


9 (Rx)  ■ . 9 (Rz)  . 9 Rx  , , 9 Rz  . , „ 

D . . = [ 3 liiL  . _--xL?kl  * 

d»k  k 3z,  9z,  " 9X 


‘j,k 


‘d»k 


“dtk 


For  the  SOR  iterations,  x£+?  and  z^*1  are 

0 » & J > ^ 

by  point  as 


calculated  point 
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£+ 1 
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(2.57) 


where  x*>k  and  z*>lc  are  the  updated  values  obtained  on  the 
current  iteration  of  Eqn.  (2.55).  The  quantity  w is  a pure 
number  in  the  range  0 < w < 2,  which  is  called  the 
relaxation  factor. 

The  application  to  the  one-dimensional  case  is 
similar.  Thus,  the  residual  error  of  Eqn.  (2.48)  is  defined 
as 


(Rs),  = 

J 


(S?5  +(— )WoS^/[2(1+  -*-S^W)]  },  , 


X ' 
w 


s K 


,5 


w 


(2.59) 


j = 2,  . . .,  M-1 


The  Taylor  series  expansion  about  S.  gives 

J 


<->r  • (rs)‘+(s‘”-s*)[  8S 


(2.60) 


o „ 3(Rs) . 

where  (Rs)  . = Rs(S.)  and  the  derivative  term,  9g  a,  is 

J J 

simply  calculated  by  differentiating  the  discretized  central 

difference  form  of  Eqn.  (2.59) » 


9[(Rs)j]  /3Sj  = 


-2, 


(2.61) 


34 


Setting  (Rs)^+1  = 0 in  (2.60)  and  solving  for  Sf+1  point  by 
J J 

point  give 

Sf+1  = S*  + AS*  (2.62) 

J J J 

where 

ASf  = -^(Rs)j’  (2.63) 

Finally,  the  N-R  method  with  the  SOR  procedure  added  into 
Eqn.  (2.62)  is  formed  as 


S*+1  = S*  + u>[(S*)*+1-S*] 
J J J J 


(2.64) 


0 4-  i 

in  which  (S*);:  is  the  updated  value  obtained  on  the 

J 

current  iteration  of  Eqn.  (2.62).  The  factor  u may  be 
specified  as  the  same  quantity  as  used  in  Eqn.  (2.57). 

This  method  can  be  implemented  by  an  alternating 
direction  such  that 


Iterations  = Odd  number  Iterations  = Even  number 


II. 5. 2 Implicit  Line  Iteration  (I.L.) 

This  method  is  to  apply  a block  iterative  method 
implemented  line  by  line.  In  order  to  do  this,  the  left 
hand  side  of  Eqn.  (2.49)  needs  to  be  assembled  into  block 
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matrix  form.  For  example,  in  Eqn.  (2.49),  let  us  define  A = 
^(a^+a^),  B = -2(b>j+b^)f  and  C = -2(c>|+c^),  also  set  M = N 
= 6 such  that  j=1,2,  ...,6;k=1,2,  ...,6.  The 
block  matrix  is  readily  formed  in  the  following  block  tri- 
diagonal  form  for  the  kth  row: 
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(2.65) 


k = 2, 


, N-1 
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The  right  hand  side  of  Eqn.  (2.65)  contains  knovm  values 
from  the  previous  iteration  " z"  and  the  boundary  condi- 
tions. Computational  algorithm  implemented  row  by  row,  and 
the  unknown  values  are  solved  by  linearlizing  their  coeffi- 
cients in  the  block  matrices.  Similarly,  for  column, 
Eqn.  (2.50)  can  also  be  formed  and  executed  column  by  column 
to  obtain  the  solutions,  such  that  the  form  is 
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a3  b3 
b3  o3 


a3  b3 
b3  °3 


A B 
B C 


a3  b3 
b3  °3 


a3  b3 
b3  °3 


A B 
B C 


a3  b3 
b3  °3 


a3  b3 
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(2.66) 
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<j  — 2 , . . • , M-1  . 

II. 5. 3 Alternating  Direction  Implicit  (ADI) 

We  combine  Eqn.  (2.65)  and  (2.66)  by  implementing  a 
first  iteration  for  line  in  the  row  direction  and  then  a 

second  iteration  for  the  j*'*1  line  in  the  column  direction 
for  a complete  cycle  of  iteration.  Such  methods  are  aptly 
designated  "alternating  direction  implicit"  (ADI)  methods. 
The  process  of  ADI  is  briefly  described  as  follows: 


II .6  Numerical  Experiments  and  Results 
The  Harris  800  computer  of  the  engineering  college  of 
the  University  of  Florida  was  employed  to  perform  the 
computations.  The  numerical  experiments  were  performed  to 
find  an  optimal  relaxation  factor  to.  These  experiments  also 
examine  the  effects  of  orthogonality  and  concentration,  and 
analyze  the  efficiency  and  accuracy  of  the  results.  Two 
experimental  problems  have  the  same  size  domain  with 
rectilinear  boundaries  and  12x12  grid  points  (i.e.,  M = N = 
12).  However,  the  distribution  of  the  grid  points  on  the 
boundaries  is  different;  for  the  first  problem,  the  grid 
spacing  is  an  arbitrary  choice;  for  the  second,  the  grid 
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spacing  is  uniform.  In  the  physical  space,  the  ranges  of  x 
and  z are  specified  by 


x : [ x1 , xM]  = [1,12] 
z : [zvzN]  = [1,12] 


(2.67) 


The  allowable  error,  RQ,  is  defined  as  the  error  of 
allowable  residual.  For  convenience,  the  finite  difference 
form  of  residual  Eqn.  (2.51)  is  rewritten  here  as 


(Rx).>k  = [LA(x)+LB(z).»u^|<i)2]J>k 


(2.68) 

(Rz)d>k  - [LB(x)+LC(z)tAu%i)2JJik 

where  j = 2,  3,  • . . , M-1 ; k = 2,  N-1 . If  every 

point  (j,k)  satisfies  the  conditions 


(Rx)j,k  — Ro 


(Rz)j,k^Ro 


then  we  find  an  approximate  solution  with  this  specified 
allowable  error  RQ.  The  number  of  iterations  and  the  amount 
of  CPU  time  required  to  reach  an  approximate  solution 
satisfying  these  conditions  at  every  grid  point  are 
determined  by  the  numerical  experiment. 

Since  the  three  methods  are  implemented  with  SOR,  it  is 
desirable  to  find  an  optimal  relaxation  factor,  in  order  to 


39 


obtain  the  highest  efficiency  for  each  method.  The  first 
problem  is  investigated  for  this  purpose.  The  results  are 
shown  in  Fig.  2.3,  which  illustrates  the  relation  between 
the  number  of  iterations  (N)  and  the  relaxation  factor 
U).  These  results  are  obtained  with  X0  = xw  = 1 and  the 
maximum  allowable  error,  RQ  = 10”^,  for  every  point  within 
the  domain.  These  results  obtained  by  the  N-R  method  in 
Fig.  2.3(a)  and  by  I.L.  method  in  Fig.  2.3(c)  clearly  show 
the  sensitivity  of  the  number  of  iterations  required  to  the 
value  of  to.  The  N-R  and  ADI  [Fig.  2.3(b)]  methods  give 
wider  range  of  relaxation  factor  than  the  I.L.  method.  The 
optimal  relaxation  factor  of  N-R,  ADI,  and  I.L.  for  this' 
experiment  is  found  as  1.6,  1.5,  and  1.45,  respectively, 
based  on  the  lowest  number  of  iterations  from  each  curve. 
Also,  this  example  problem  is  studied  for  the  effects  of 
smoothness,  orthogonality,  and  concentration, 
respectively.  The  weight  function  is  specified  as 

W(x,z)  = (x+z) [24-(x+z) ] . (2.69) 

Figure  2.4(a)shows  the  result  of  minimizing  the  smoothness 

functional  only  (I  = I , x = ^ — = 0,  x = ^ — =0);  it 
J n s’  Ao  - , w ~ , 

Ao  Aw 

shows  the  grid  points  distributed  smoothly.  Figure  2.4(b) 
gives  a result  of  minimizing  the  functionals  of  smoothness 
and  orthogonality  (XQ  =8,  Xw  = 0).  Clearly,  the 
orthogonality  near  the  boundary  is  more  strongly  emphasized; 
however,  the  smoothness  is  not  as  good  as  in  Fig.  2.4(a). 
Figure  2.4(c)  shows  a result  of  minimizing  the  functionals 
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Fig.  2.3 


Relaxation  factor  ( u>)  versus  number  of  iterations 

(N) 
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Fig.  2.4  Grid  generations  with  unequal  spacing  on  boundary 
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of  smoothness  and  concentration  (XQ  =0,  Xw  = 8).  Since  the 
maximum  values  of  W(x,z)  are  on  the  line  x+z  = 12  [see  the 
dashed  line  in  Fig.  2.4(c)],  it  emphasizes  the  concentration 
more  strongly  than  the  region  which  is  away  from  this 
line.  Also,  the  grid  in  Fig.  2.4(c)  is  skewed  in  comparison 
with  the  grid  in  Fig.  2.4(a);  however,  if  the  adaptive 
boundary  grid  were  applied  to  this  boundary,  the  skewness 
would  be  reduced.  Thus,  the  different  control  parameters 
result  in  different  properties  of  the  grid  distributions. 
More  detailed  descriptions  concerning  the  interaction  of 
these  parameters  and  changes  can  be  found  in  references  [6] 
and  [19] . 

The  second  example  problem,  grid  generation  with 
uniform  grid  spacing,  examines  the  influence  of 
concentration  on  the  grid  network.  The  weight  function  is 
specified  as 

W(x,z)  = (x+z) [26-(x+z) ] (2.70) 

The  maximum  values  of  W(x,z)  are  on  the  line  x+z  = 13  and 
W(x,z)  is  decreased  away  from  this  line.  Large  values  of 
W(x,z)  cause  higher  concentration,  smaller  cells,  and  small 
values  cause  lower  concentration.  Figure  2.5(a),  (b) , and 
(c)  show  the  results  of  minimizing  the  functionals  of  the 
smoothness  with  their  concentrations  Xw  = i , io,  and  100, 
respectively.  The  grid  generated  in  Fig.  2.5(b)  is  clearly 
different  from  Fig.  2.5(a),  because  larger  values  of  the 
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(a)  A/x^  = o,  A/x;  = 1 


(c)  a/x^  = o,  x/x^  = 100 

Fig.  2.5  Grid  network  with  different  A and  uniform  grid  on 
boundary 
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specified  xw  place  more  emphasis  on  concentration.  However, 
the  difference  between  Fig.  2.5(b)  and  2.5(c)  is  not  as 
pronounced.  This  indicates  that  the  variation  of  the 
concentration  functional  converges  to  some  constant  value 
(i.e.,  WJ“2  = const.)  when  xw  >>  1.  In  Fig.  2.5(c),  the 
ratio  of  maximum  value  of  W(x,z),  Wmax(x,z)  = 169,  to 
minimum  value,  Wmin(x,z)  = W(2,2)  = W(11 ,11)  = 88,  is 
169/88.  The  corresponding  ratio  of  minimum  to  maximum 
values  of  cell  area,  Amin/^max  = (83/169)^^  a 0.72,  is 
close  to  1 that  is  why  the  grid  sizes  near  and  far  from  the 
line  x+z  = 13  may  not  be  clearly  apparent  in  Fig.  2.5- 

In  order  to  analyze  the  efficiency  and  accuracy  of  the 
results,  the  average  error,  E,  is  defined  as 

E = [ e(x-x  ) 1+  E ( z — z ) . ]/2[(M-2)(N-2)]  (2.71) 

O J , K O J , K. 

where  j = 2,  3,  - - . , M— 1 ; k = 2,  3,  - - . , N-1 ; the 
notation  z represents  the  summation  on  both  j and  k;  and 
(x0,zQ)jjk  are  the  convergent  solutions  found  with  the 
allowable  error,  RQ  = 10“^,  while  the  approximate  solutions 
(x,z)j>lc  correspond  to  some  (considered  convergent)  larger 
allowable  RQ.  The  first  problem  with  XQ  = = 1,  xw  = 

x/ = 1 is  used  for  this  experimental  analysis.  Figure  2.6 
illustrates  the  efficiency  analysis  for  those  methods,  which 
have  the  property  that  CPU  time  is  linearly  related  to  RQ; 
the  slope  of  the  line  of  the  Jacobi  Method  (without  update 
values  and  SOR  procedure  in  N-R  method)  is  highest 
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Ro 


Fig.  2.6  Efficiency  analysis 
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and  N-R  with  SOR  results  in  the  lowest  value  of  CPU  time  and 
slope  in  comparison  with  the  others.  Figure  2.7  illustrates 
the  relations  between  average  error  E and  R0,  it  shows 
clearly  that  the  larger  differences  in  the  values  of  E for 
the  different  methods  occur  at  RQ  = 1 0~ ^ and  that  the 
differences  are  smaller  up  to  RQ  = 10”^;  N-R  method  with  SOR 
results  in  the  lowest  E over  all  RQ  up  to  10-^.  However, 
usually,  for  most  of  the  large  number  of  grid  generation, 
the  allowable  error  RQ  up  to  10-4  is  good  enough,  since  the 
grid  generation  does  not  need  to  reach  too  highly  accurate 
grid  solution;  also,  a large  CPU  time  would  be  expended  to 
obtain  a greater  accuracy.  In  the  view  of  efficiency  and 
accuracy,  we  conclude  that  the  N-R  method  with  SOR  is  the 
most  suitable  of  the  methods  examined  to  be  used  for  the 
adaptive  grid  generation. 


CHAPTER  III 

NAVIER-STOKES  EQUATIONS 


III.1  Review  of  Fundamental  Principles 


The  purpose  of  this  section  is  to  review  the  set  of 
governing  equations  resulting  from  the  physical  laws  of 
conservation  of  mass,  momentum,  and  energy  for  three- 
dimensional,  compressible  fluid  flow. 

The  continuity  equation,  a consequence  of  the 
conservation  of  mass,  is 


where  p = p(x,y,z,t)  is  density,  u = u(x,y,z,t),  v = 
v(x,y,z,t),  and  w = w(x,y,z,t)  are  the  Cartesian  velocity 
components.  The  notations  3X  = 3/ax,  . . .,  etc.,  are 
partial  derivatives. 

The  Navier-Stokes ' equations,  a consequence  of  the 
conservation  of  momentum,  are 


(3.1) 
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in  which  = d/dt  is  a substantial  derivative  which 

consists  of  a local  and  convective  contribution,  p is  the 
pressure,  and  the  stress  components  for  Newtonian  fluid  are 
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where  X = - — \i  (no  bulk  viscosity)  and  u are  viscosity 
coefficients.  The  subscripts  x,  y,  and  z in  the  right  hand 
side  of  Eqn.  (3*3)  denote  the  partial  differentiations  with 
respect  to  x,  y,  and  z,  respectively. 

The  energy  equation,  a consequence  of  the  conservation 
of  energy,  can  be  written  in  the  form 
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(3.4) 


where  e is  the  total  energy  per  unit  volume,  given  by 
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e = p[eI  + ^-(u2+v2+w2)]  , (3.5) 

and  ej  is  the  internal  energy  per  unit  mass;  < is  the 
coefficient  of  thermal  conductivity,  T is  the  temperature 
and  Tx,  Ty,  and  Tz  are  the  partial  derivatives  of  T with 
respect  to  x,  y,  and  z,  respectively. 

For  the  perfect  gas,  the  thermal  and  caloric  equations 
of  state  can  be  defined  as 

p = pRT , eI  = cyT,  (3.6) 

respectively,  where  R is  the  gas  constant,  and  cv  is  the 
specific  heat  at  constant  volume.  Also,  the  following 
relationships  exist 


h = 
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R 

Y-1 


YR 

Y-1 


(3.7) 


where  h is  the  enthalpy,  Cp  is  the  specific  heat  at  constant 
pressure  and  y is  the  ratio  of  specific  heats.  For  air  at 
standard  conditions,  y = 1.4.  Thus,  the  perfect  gas 
equation  of  state  can  also  be  written  as 


P = (Y-1)[e  - -^-p  (u^+v^+w^)  ] 


(3.8) 


The  coefficients  of  viscosity  and  thermal  conductivity 
have  been  related  to  the  thermodynamic  variables  using 
kinetic  theory.  It  is  possible  to  use  an  interpolation 


formula  based  on  D.M.  Sutherland's  theory  of  viscosity  [21] 


(3.9) 


where  denotes  the  viscosity  at  the  reference  teraprature 
T,,,,  and  S is  the  Sutherland  constant,  which  for  air  assumes 
the  values  S = 193. 6°R  or  110°K.  The  Prandtl  number, 

c u 

Pr  = , (3.10 

is  often  used  to  determine  the  coefficient  of  thermal 
conductivity  < once  u is  known.  This  is  possible  because 
the  ratio  (c^/Pr)  which  appears  in  the  expression 


is  approximately  constant  for  most  gases. 

The  derivatives  of  ej  in  the  caloric  equation  (3.6) 
with  respect  to  x,  y,  and  z at  cv  = constant  can  be  written 
as  the  vector  form, 
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K = — 2.  p 


(3.11) 


I 


After  substituting  this  equation  and  Eqn.  (3.11)  into  Eqn. 
(3-4)  and  shifting  the  pressure  terms  to  the  left  hand  side, 
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the  energy  equation  can  be  rewritten  as  the  form 
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(3.13) 


III. 2 Conservation-Law  Form  of  Equations 
A partial  differential  equation  (PDE)  can  be  put  in 
conservation-law  form,  if  the  coefficients  of  the  derivative 
terms  are  either  constant  or  if  variable,  their  derivatives 
disappear  in  the  equation.  In  either  case,  the  coefficients 
can  be  brought  inside  the  derivative.  Normally,  physical 
conservation  statements  are  represented  by  PDEs;  this  means 
that  the  divergence  of  a physical  quantity  can  be  identified 
in  the  equation  [22].  Using  an  appropriate  set  of  equations 
in  conservation-law  form  represented  by  the  PDEs  of  interest 
which  all  have  their  basis  in  physical  laws  such  as  the 
conservation  of  mass,  momentum,  and  energy  shows  the  nature 
of  the  physical  conservation  laws  involved  and  allows 
overall  integral  properties  to  be  maintained  identically  in 
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the  finite-difference  system.  Such  a system  of  equations  is 
now  in  common  use  in  calculations  of  very  high  gradient  for 
a physical  phenomenon,  regardless  of  the  finite-difference 
methods  used,  since  the  exact  planar  high  gradient  property 
such  as  shock  speed  will  be  produced  by  any  stable  method 
[23,24]. 

The  continuity  equation  (3*1)  and  energy  equation 
(3.12)  are  already  in  conservation-law  forms.  The 
conservation-law  form  of  the  momentum  equation  can  be 
readily  formed  by  adding  the  continuing  equation  (3.1)  into 
Eqn.  (3.2).  The  result  is 
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Also,  the  compressible  Navier-Stokes  equations  in  Cartesian 
coordinates  can  be  formed  by  assembling  Eqn.  (3.1),  (3.14), 
and  (3-12)  into  the  compacted  form 


V + 


a f + 

X X 


9 f + 

y y 


9 f 
z z 


0 


(3.15) 
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where  q,  fx>  f^,  and  f^  are  the  five-component  vectors 

IT1 

q = { p pu  pv  pw  e } 
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The  governing  equations  can  be  made  dimensionless.  The 
viscosity  the  density  Pao,  the  temperature  T^,  the  sound 
speed  aOT,  and  a characteristic  length  D are  chosen  as 
reference  quantities,  where  the  subscript  ® denotes  the  free 
stream  condition.  Let 


(3.16) 


in  which  the  superscript  * denotes  the  nondimensional 
properties.  Substituting  Eqn.  (3.16)  into  Eqn.  (3.15)  ana 
finally  dropping  the  superscript  * lead  to  the  following 
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dimensionless  form  of  the  equations 


3tq  + 9x(E"Ev)  + 8y(F“Fv)  + 3z(G“Gv)  = 0 (3.17) 

where  q,  E,  F,  and  G are  five-component  vectors, 
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and  E , Fv,  and  Gv  are  the  viscous  flux  terms, 


E = Re"V  x 
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(3.19) 


The  Reynolds  number  Re  is  defined  as 
= p a D/ y 

CO  00  00 


Re 


(3.20) 
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III. 3 Transformed  Governing  Equations 
As  mentioned  previously,  it  is  possible  to  develop  the 
well  known  "boundary-fitted  curvilinear  coordinate  systems" 
by  coordinate  transformation  so  that  unequal  spacing  in 
physical  domain  can  be  mapped  onto  equal  spacing  in  the 
computational  domain.  The  governing  equations  need  to  be 
transformed  in  the  same  manner.  With  this  procedure,  the 
solution  of  a partial  differential  system  is  calculated  in  a 
fixed  parallelepiped  field  with  a cubic  grid  in  the 
computational  domain.  Also,  the  significant  aspect  of  the 
transformation  offered  in  Chapter  II  is  that  grid  points  can 
be  adaptively  clustered  in  regions  that  have  rapid  changes 
in  the  flow  field  gradients. 

Let  us  now  consider  a general  transformation  of  the 
three-dimensional  form 

5 = 5(x,y,z,t) 

n = n(x,y,z, t)  (3.21 ) 

S = ?(x,y,z, t) 

T = t 

where  x is  the  computational  time.  It  can  then  be  used  to 
transform  the  governing  equations  from  the  physical  domain 
(x,y,z,t)  to  the  computation  domain  (S,ti,c,t).  The  three- 
dimensional  version  of  Eqn.  (2.26)  is 
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(3.22) 


therefore,  the  metrics  g , n , . . .,  etc.  are  found  by  the 

X X 


cofactor  method, 
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_ *1 

where  J , the  3-D  version  of  the  Jacobian  determinant 
defined  by  Eqn.  (2.15),  is  given  by 


-1  _ 3(x,y,z) 

3(5, n,0 


X X 

5 n c 


y€  yn  yc 
Z5  Zn  Zc 


(3.24) 


= x5(ynzrycZn)-xn(y5ZrycZC)+X?(y5Zn-ynZ5)  (3‘25) 


Thus  J 


represents  the  cells  volume  in  3-D  case. 


The 
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metrics  can  be  readily  obtained  after  finding  the  inverse 
transformation, 


x — x(c,n,?,x) 

y = y(5,nfc,T)  (3.26) 

z = z( £ , n , c , x) 
t = t 


If  the  (5,n,c,t)  system  is  obtained  by  solving  the  governing 
equations  of  grid  generation  (i.e.,  Euler's  equations,  as 
described  in  Chapter  II  Adaptive  Grid  Generation),  we  may 
solve  any  set  of  PDE  on  this  physical  coordinate  system  by 
solving  the  transformed  equations  in  the  computational 
domain  [25].  Also,  the  transformed  computational 
coordinates  of  a given  moving  grid  point  in  the  physical 
space  are  fixed;  thus, 


dg  _ dri_  _ d£ 
dx  " dx  ~ dx 


(3.27) 


gives 


where  the  subscripts  t and  x denote  partial  differentiations 
with  respect  to  physical  time  t and  computational  time  t, 
respectively. 
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We  apply  this  generalized  transformation  to  Eqn.  (3.17) 
and  obtain  the  following  transformed  equation 

qT  * «tq5  + "tq„  + ?tq? 

+ SX(E-Ev)?  * * Cx(E-Ev)5 

(3.29) 

* 5x(F-Fv)5  . nx(F-Fv)n  . CX(F-Fv)t 

+ + "x(G‘Gv)„  + 'x(G*Vc  = 0 

This  equation  may  be  written  with  the  metric  coefficients 
inside  the  differentiation,  thus  producing  a source  term. 

The  resulting  form  is  called  the  "weak  conservation  law 
form" 


q t * 3C[5tq  * 5X(E-EV> 

+ 5n[ntq  * nx(E_Ev) 

+ 3cI«tq  + «x(E_V 


- «y(F-Fv)  ♦ ?Z(0-GV)] 

+ ly(F-Fv)  + nz(G-Gv)] 

+ Cy(F-Fy)  + cz(G-Gv))  + S' 


(3.30) 


0 


where 


S' 


[(5t}5  + (nt)T1  + 


+ (ny)n 


■[(Cz)e  + ^zK  + 
z t,  z n 


(?t)c]q  - [(Sx)?  - Ux)„  * (CX){](E-EV) 
(Cy);](F-Fv) 


It  has  been  shown  by  Viviand  [26]  and  Vinokur  [27]  that  Eqn. 
(3.30)  can  be  put  back  into  "strong  conservation-law 
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form."  In  order  to  do  this,  Eqn.  (3.29)  is  first  divided  by 
the  Jacobian  and  is  then  rearranged  into  conservation-law 
form  by  adding  and  subtracting  like  terms.  The  following 
equation  results 

aT(J-1q)  + 35{J-1[5tq  + ?X(E-EV)  + Sy(F-Fv)  + ?2(G-Gv)]} 

+ [ntq  + VE-Ev)  + ny(F-Fy)  + nz(G-Gy)]} 

+ 9^J‘1[?tq  + Cx(E_Ev)  + Cy(F“Fv}  + ?2(G_Gv)]^  + S = 0 

(3.31) 


where  the  source  terms, 

- [(J_15x)5  ♦ - (J-1Cx)c](E-Ev) 

- l(J‘15y)5  + (J‘\>n  + (J_1Cy)5](F-Fv) 

- [(J’15z)|  ♦ (J't2)n  - (J'1CZ)?](G-GV)  , 

are  all  equal  to  zero  and  can  be  dropped.  This  can  be 
verified  by  substituting  the  metrics  given  by  Eqn.  (3.23) 
into  this  form.  Therefore,  Eqn.  (3.31)  leads  to 

9Tq  + 3 £ (E-Ev)  + 3n(F-Fv)  + 3?(G-GV)  = 0 


(3.32) 
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where  q,  E,  F,  and  G are 


r-1 


p 

PU 

pV 

pu 

puU+£xP 

PuV+nxp 

PV 

“ -1  I 

> > E = J S 

pvU+5yP 

1 

►"D 

II 

< U-* 

pvV+Tiyp 

PW 

pwU+5zP 

pwV+nzP 

e 

\ j 

(e+p)U-5tP 
^ - 

(e+p)V-ntp 

«.  j 

>> 


~ -1 
G = J < 


pW 

puW+?xp 

pvW+Cyp 

pwW+?zp 

(e+p)W-etp 


(3.33) 


in  which  U,  V,  and  W are  contravariant  velocities  defined  as 
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(3.34) 


and  the  transformed  viscous  flux  terras  Ev,  Fy,  and  Gy  are 
given  by 


O » **3  > M > 
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CHAPTER  IV 

THIN-LAYER  NAVIER-STOKES  EQUATIONS 
AND  COMPUTER  CODE 

IV. 1 Thin-Layer  Navier-Stokes  Equations 
In  the  thin-layer  approximation  to  the  Navier-Stokes 
equations,  the  viscous  terms  containing  derivatives  in  the 
directions  parallel  to  the  body  surface  are  neglected  in  the 
unsteady  Navier-Stokes  equations,  but  all  other  terms  in  the 
three  momentum  equations  are  retained.  One  advantage  of 
retaining  those  terms  which  are  normally  neglected  in 
boundary-layer  theory  in  the  momentum  equations  is  that 
separated  and  reverse  flow  regions  can  be  readily  com- 
puted. Also,  flows  which  contain  a high  normal  pressure 
gradient,  such  as  the  flow  fields  where  the  boundary-layer 
equations  are  not  applicable,  can  be  calculated  [22]. 

Steger  [10]  was  the  first  to  consider  a two-dimensional 
thin-layer  model  for  simulation  of  flow  about  arbitrary 
geometries  with  application  to  airfoils  problem.  Later, 
Baldwin  and  Lomax  [28]  evolved  the  thin-layer  approximation 
from  a detailed  investigation  of  very  high  Reynolds  number 
computation  from  the  full  Navier-Stokes  equation.  In  these 
computations,  a substantial  fraction  of  the  available 
computer  storage  and  time  is  expended  in  resolving  the 
normal  gradients  in  the  boundary  layer  since  a 
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highly  concentrated  grid  is  required.  As  a result,  the 
gradients  parallel  to  the  body  surface  are  usually  not 
resolved  in  an  adequate  manner  even  though  the  corresponding 
viscous  terms  are  retained  in  the  computations.  Hence,  for 
many  Navier-Stokes  computations,  it  makes  sense  to  drop 
those  terms  provided  that  they  are  reasonably  small.  It 
should  be  emphasized  that  the  thin-layer  approximation  is 
valid  only  for  high  Reynolds  number  flow. 

For  the  application  of  the  thin-layer  approximation,  it 
becomes  necessary  to  map  the  physical  body  surface  onto  a 
transformed  coordinate  surface  (i.e.,  5 = const.).  Thus, 
the  transformed  governing  equation  (3*32)  can  be  simplified 
by  dropping  all  viscous  derivatives  in  the  direction 
parallel  to  the  body  surface  (c  = const.),  and  we  obtain  the 
thin-layer  Navier-Stokes  equations, 

A A A A yj  A 

3q+3E+3F+3G=Re  9S  (4.1) 

t 5 n 5 4 

AAA  A 

where  q,  E,  F,  and  G were  defined  by  Eqn.  (3. 33)  while  the 

A 

viscous  terms  S are  readily  obtained  by  retaining  the 

A 

viscous  flux  terms  in  the  c-direction  in  Gv  of  the  Eqn. 
(3.35). 
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(4.2) 


IV. 2 Surface  Boundary  Conditions 


For  the  boundary  conditions  along  the  surface 
?(x,y,z,t)  = const.,  the  velocities  u,  v,  and  w are  obtained 
from  Eqn.  (3.34),  such  that 


In  inviscid  flow,  W can  be  defined  as  W = 0,  while  U and  V 
are  unknown  values  which  can  be  found  by  extrapolation  from 
interior  points.  The  pressure  along  the  body  surface  can  be 
obtained  from  a normal  momentum  relation  by  combining  those 
three  transformed  momentum  equations. 


(4.3) 
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(4.4) 


where  Pn  is  the  normal  pressure  gradient  at  the  body  surface 
and  P^,  P^,  and  P?  are  the  pressure  gradients  along  the 
directions  of  £,  n,  c,  respectively. 

For  viscous  flows,  the  equations  (4.3)  and  (4.4)  are 
used  with  U = V = W = 0.  All  of  the  preceding  boundary 
conditions  are  valid  for  steady  or  unsteady  flow. 


The  effects  of  turbulence  can  be  simulated  in  terms  of 
the  eddy  viscosity  coefficient  y^..  The  effective 
coefficient  of  viscosity  yg  is,  therefore,  defined  by 


where  y is  the  dynamic  laminar  viscosity.  Recently,  the 
two-layer  eddy  viscosity  closure  model  has  been  extensively 
used  in  solving  practical  turbulent  flow  problems.  This 
model  consists  of  inner  and  outer  regions,  consequently,  the 
distribution  of  the  two-layer  effective  eddy  viscosity 
across  the  boundary  layer  can  be  written  as 


IV. 3 Turbulence  Model 


ye  = y + yt 


(4.5) 
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(4.6) 


(we)o  = y + (yt}o 


in  which  subscripts  i and  o represent  the  inner  and  outer 
regions,  respectively.  Similarly,  the  effective  heat  flux 
term  ue/(Pr)e  can  be  defined  as 


where  (Pr)e  and  (Pr)-j-  are  viewed  as  effective  and  turbulent 
Prandtl  numbers,  respectively.  Equations  (4.6)  and  (4.7) 
are  applied  to  the  transformed  thin-layer  Navier-Stokes 
equation;  thus,  the  y and  Pr  in  Eqn.  (4.2)  are  replaced  by 
ue  and  (Pr)e,  respectively. 

The  turbulence  closure  model  programmed  in  the  thin- 
layer  Navier-Stokes  code  is  an  extension  of  the  Cebci’s 
algebraic  eddy  viscosity  model  for  avoiding  the  finding  of 
the  boundary  layer  edge  [28].  This  model  in  terms  of  the 
two-layer  algebraic  eddy  viscosity  is  given  by 


e 


(4.7) 


(^t^i  = p{0.4z[1-exp(-z+/A+)  ] }2  iwi  , z _<_  zc 


(4.8) 


(y.)  = 0.0168PC  F 

to  cp 


F F 

cp  wake  kleb 


(z)  , z zc 


where  subscripts  i and  o represent  the  inner  and  outer 
regions,  respectively.  In  Eqn.  (3.16),  z and  p are  the 
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dimensionless  quantities,  the  normal  distance  from  the  wall 
and  the  density,  respectively;  zc  is  the  smallest  crossover 
value  of  z at  which  values  from  inner  and  outer  formulas  are 
equal.  The  constant  A+  denotes  an  empirical  constant  equal 
to  26.  The  law-of-the-wall  coordinate  z+,  vorticity  mu  , 
Klebanoff  intermittency  factor  Fkleb(z),  wake  factor  Fwake, 
and  additional  constant  were  defined  in  reference 
[28] . This  model  was  applied  to  obtain  the  comparisons 
which  were  made  with  experiment  for  an  incident  shock  on  a 
flat  plate,  separated  flow  over  a compression  corner,  and 
transonic  flow  over  an  airfoil.  The  results  show  that 
separation  and  reattachment  points  from  numerical  Navier- 
Stokes  solutions  agree  with  experiment  within  one  boundary- 
layer  thickness.  Also,  Pr  and  (Pr)t  are  defined  by  Pr  = 
(Pr)^  = 0.9  as  reported  in  [28]. 

IV. 4 Computer  Code 

As  described  in  the  Introduction,  Nietubicz,  Pulliam, 
and  Steger  [14]  developed  a numerical  procedure  for  an 
axisymmetric  case  form  of  the  thin-layer  Navier-Stokes 
equations.  This  axisymmetric  thin-layer  code  with 
cylindrical  coordinate  system  which  has  readily  been 
applied  to  transonic  projectile  aerodynamics  problems  at 
the  U.S.  Army  Ballistic  Research  Laboratory  is  more 
efficient  than  the  thin-layer  Navier-Stokes  Code  with 
Cartesian  coordinate  system.  A sketch  of  a typical 
axisymmetric  body  is  shown  in  Fig.  4.1(a).  In  order  to 
determine  the  circumferential  variations  of  typical  flow  and 
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Z t 


(a)  Projectile  type  body 


(b)  x = const,  plane 


R 


(c)  <t>  = const,  plane 


Fig.  4.1  Axisymraetric  body  and  coordinate  system 
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geometric  parameters,  we  first  establish  correspondence 
between  the  Cartesian  coordinates  (x,y,z),  the  cylindrical 
coordinates  (x,sz$,R)  in  the  physical  space,  and  the 
transformed  variables  (C,n,?).  For  the  cylindrical 
coordinates,  it  may  be  useful  to  have  the  reference  line  for 
<t>  rotate  relative  to  the  (x,y,z)  system,  so  that  £ = *5(x), 
where  t is  the  time  in  (£,n,s)  system.  For  example,  the 
cylindrical  coordinate  grid  system  may  be  chosen  to  rotate 
relative  to  a spinning  body.  From  the  views  shown  in  Fig. 
4*1  > the  relationship  between  the  coordinate  systems  are 
observed  to  be 

<b  = on 

x = x(  5, 5,  t) 

y = RU,s,T)sin?$  (4.9) 

z = R(  5 , ? , T)cosg$ 
t = t 

where  c is  an  arbitrary  constant  which  can  be  set  equal  to 
one,  and  £ is  defined  as  6 = tf(x).  Note  that  x and  R are 
functions  of  £,  £,  and  time  x only.  Thus,  the  metric  terms 
are  readily  evaluated  from  Eqn.  (4.9),  such  that 


'XS  Xn  X?  xx  ' 

' X5  0 xc 

yc  yn  yc  yT 

= 

R^sin^  s^Rcosjri  R^sin^  RTsin$+^TRcosj$ 

.ZS  ZX|  ZS  ZT  . 

R^cosd  -<z$nRsin$  R^cos$  RTcos^-^TRsin6 

(4.10) 
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Now  substituting  Eqn.  (4.10)  into  Eqn.  (3.23)  and  (3.28) 
yields 


5x  «y  5z  5t 
nx  ny  nz  nt 
?x  5y  ?z  ?t 


-1 


jzS^RR  ^ s^Rx^sintf 
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-sin^/  ( JR^zS ^ ) 

-iTnun) 

-SZ^RR^  j^Rx^sin^ 

i>  Rx-cosi i 
"n  s 

«iI1R(xTRrRtxs) 

(4.11) 


in  which  the  Jacobian  is  given  by 


J-1  = R(z5n(x5Rc-x?R5)  (4.12) 

a 

After  the  term  is  evaluated  in  Eqn.  (4.1),  we  can 
eliminate  the  cylindrical  coordinates  altogether  by  letting 
f>  = 0°,  hence  R = z (see  Fig.  4.1b),  sin^  = 0 and  cos(Z$  = 

1.  The  term  is  simply  the  scaling  constant,  jZ$n  = 1 , and 
the  Cartesian  velocities  defined  in  Eqn.  (4.3)  become 
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It  should  be  emphasized,  for  an  axisymmetric  system,  the 
state  variables  and  the  contravariant  velocities,  U,  V,  and 
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¥ are  required  to  be  invariant  in  the  ^direction.  The 
metrics  terms  n * and  ^ are  zero  for  steady  non- 
spinning projectile  flows  and  also  are  invariant  in  the  n- 
direction  for  unsteady  flow  caused  by  body  motion  such  as 
axial  acceleration  or  axisymmetric  spinning. 

The  resulting  expressions  for  u,  v,  w,  and  the  metrics, 
Eqn.  (4.11),  have  only  sinjzi  and  cos^  variation  in  the  n- 
direction.  Consequently,  Eqn.  (4. 11)  and  (4.13)  can  now  be 
used  in  Eqn.  (4.1)  to  reduce  the  n derivative  of  the  flux 
vector  F [Eqn.  (3-33)],  such  that 
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Transformation  of  the  n derivative  results  in  the  following 
term: 
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The  resulting  thin-layer  n-invariant  equations  are  written  as 
3 q + 3 _E  + 3 G + H = Re_13rS  (4.16) 


with  the  metrics  and  Cartesian  velocities  in  q,  E,  G,  and  S 
of  Eqn.  (4.1)  evaluated  at  <t>  = 0°,  and 


H 


J"V 


0 

0 

i PV[R5(U-5t)+R5(W-Ct)]  > 


-PVR^n(V-nt)-p/(R^n) 


0 


(4.17) 


For  surface  boundary  condition,  the  body  surface  pressure 
defined  in  Eqn.  (4.4)  is  simplified  as 


P (c2+<;2)2 

nv  x z 


(5  c +5  C )Pr  + U2+S2)Pr 
v x x z z'  5 x z'  ? 


= p [8T?t+u3T4x  + v(  J^T1RXg(6T  )+w3t  ( Jj^Rx^  ) ] 


- PU(^xu5  + Czw5)+PV[Jx5R2i6j(V-nt)]  (4.18) 


where  Pn  is  the  normal  pressure  gradient  at  the  body  surface. 

For  an  axisymmetric  body  spinning  with  angular  velocity 
w,  we  have  the  choice  of  using  a (S,71,^)  curvilinear 
coordinate  grid  in  the  physical  space,  so  that  body  spins 
relative  to  the  coordinate  grid,  or  of  having  the  grid  system 
rotate  with  -u  relative  to  body.  In  the  first  case,  the  no 
slip  condition  would  be 
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U = W = 0 and  V = u (4.19(a)) 

In  the  second  case,  which  is  used  in  this  study,  the  no  slip 

boundary  conditions  is  enforced  by  setting 

U = V = W = 0 (4.19(b)) 

where  U,  V,  and  W are  contravarient  components  of  the 
relative  velocity,  relative  to  the  body.  Therefore,  for  the 
spinning  body  the  usual  no  slip  condition,  Eqn.  (4.19(b)),  is 
imposed  with  the  time  metrics  set  at  i>  = 0°  (y=0)  and 
consequently  Cy  = nz  = ?y  = 0 in  Eqn.  (4.11).  Thus,  Eqn. 
(3.28)  is  simplified  as 

5t  - .-y,VZT52  ’ -“<ZSy'y52  ’ 0 

nt  ■ -yTny-zTnz  " -“(zuy-yn2)  = -v/in  (4.20) 

st  - -ySy-zSz  - -“(^y-y^J  = 0 

Once  the  contravariant  velocities  are  specified  at  the  body, 
the  Cartesian  velocities  Eqn.  (4.13)  with  (6  = 0°  become 
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The  finite-difference  scheme  employed  here  is  the 
implicit  approximate  factorization  algorithm  used  in  the 
delta  form  as  analyzed  by  Warming  and  Beam  [12].  An  implicit 
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method  was  selected  to  avoid  restrictive  stability  conditions 
which  occur  when  small  grid  spacing  is  used.  The  solution  of 
the  two-dimensional  system  of  equations  (4*16)  is  implemented 
by  an  approximate  factorization  of  the  equations  into  two 
one-dimensional-like  systems  of  equations.  This  procedure 
has  been  employed  by  Steger  [10],  Baldwin  and  Lomax  [28], 
Pulliam  and  Steger  [11],  and  Nietubicz  [15].  Central- 
difference  operators  are  used  and  the  algorithm  produces 
block  tri-diagonal  systems  for  each  space  coordinate.  The 
stability  and  accuracy  of  the  numerical  algorithm  for 
compressible  Navier-Stokes  equations  are  described  in  detail 
by  Warming  and  Beam.  The  implicit  approximation 
factorization  algorithm  was  applied  to  the  thin-layer 
approximation  with  transformed  governing  equations  about 
arbitrary  2-D  geometries  by  Steger.  Later,  implicit  finite- 
difference  simulations  of  three-dimensional  compressible  flow 
were  developed  by  Pulliam  and  Steger. 

Using  the  Beam  and  Warming  implicit  algorithm  applied  to 
(4.16)  results  in  the  following  finite  difference  equations 

[I+hd^^-ejJ-1  A^V^J]  [l+h6^Cn-eiJ_1  A?V?J-hRe_1  6^J_1MnJ]  Aqn 

= -At[6cEn+6?Gn-Re_1 65Sn]-AtHn-£EJ_1 [(A5V5)2+(A?V?)2]Jqn 

(4.22) 

where  the  6's,  A's,  and  V's  denote  central,  forward,  and 
backward  difference  operators,  respectively,  e.g., 
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A ^ J = J(5  + A5,s)  - J(§,?).  Indices  denoting  spatial  location 
are  suppressed  for  convenience.  The  time  step  h is  defined 
as  h = At.  The  jacobian  matrics  A = 9E/9q,  C = 9G/3q  along 
with  the  coefficient  matrix  M obtained  from  the  local  time 
linearization  of  S are  described  in  detail  by  Steger  [10]. 

The  coefficients  of  the  explicit  and  implicit  "smoothing" 
terms  are  eE  and  £j,  respectively;  these  are  required  to  damp 
high  frequency  oscillations  in  the  solution.  A typical  range 
recommended  by  Nietubicz  et  al.  [14]  is  eE  = (1  to  5)At  with 

A 

eI  = 2eE.  We  define  Aqn  as 


(4.23) 


where  qn  = q(nAt)  = q(t)  is  defined  by  letting  the  temporal 
variable  t be  discretized  as  t = nAt. 


The  factored  scheme  (4.22)  for  the  system  of 
axisymmetric  thin-layer  Navier-Stokes  Equations  is, 
therefore,  implemented  as 


[ I+h^A11-^  J-”1  A^  j]  Aq*  = r.h.s.  [Eqn.  (4.22)  ] (4.24(a)) 

[I+h« _Cn-eTJ_1A  V J-hRe-1 6rJ_1anJ] Aqn  = Aq*  (4.24(b)) 


(4.24(c)) 


where  Aq*  are  intermediate  temporal  differences,  and  r.h.s. 
denotes  right-hand-side. 
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Five  subprograms,  INITIA,  EIGEN,  STEP,  MAP,  and  OUTPUT, 
are  called  in  order  in  the  MAIN  program  of  this  code.  They 
are  described  in  sequence  as  follows: 

SUBROUTINE  INITIA 

This  subroutine  computes  initial  conditions  and  their 
relative  input  parameters.  A few  input  parameters  are  read 
in  this  code;  the  detailed  descriptions  are  given  in  Appendix 
A.  The  following  set  of  constants  are  computed  before 
generating  grid  and  metrics, 

RE=Re/1M,  0MEGA=2ttw/60,  GAMMA=y=1  . 4 , SPIN=(PD0V)1M, 
GAMI=y-1,  GD=y(  Y-1 ) , HD=^At,  JM=JMAX-1 , LM=LMAX-1 , 

ND2= ( KMAX ) ( LMAX ) , DX1 =DY1 =DZ1 =1 .0,  HDX=HD/DX1 , 

HDY=HD/DY1  , HDZ=HD/DZ1  , CS=C0S ( ott/I 80 ) , and 
SS=sin(  out/  180), 

where  1M  is  the  Mach  number,  PDOV  is  a factor  transferring 
Mach  number  to  spin  number,  JMAX  is  the  maximum  number  of 
grid  points,  along  the  ^-direction,  and  a is  the  angle  of 
attack  which  must  be  set  equal  to  zero.  The  calculations  of 
grid  and  metrics  are  implemented  by  calling  the  "SUBROUTINE 
GRID"  and  "SUBROUTINE  METRIC,"  respectively.  We  then  load 

a 

the  vector  qro  that  characterizes  the  free  stream  flow 
properties  at  a = o.  Thus, 

Qoo  » { P , pu,  pv,  pw,  e}  ^ 

= {1  .0,1Mcosa,Q,1Msina,--^--  ^ + -1  1M2}£  ^ , 


• J 


l = 1 , . 


LMAX,  j=1,  . . 


JMAX 


(4.25) 
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where  the  total  flow  energy  per  unit  volume,  eOT,  can  be 
derived  from  Eqn.  (3*9), 


p.  = (Y-Ute.  - \ ml] 


pma„, 

OO  CO 


Y 


1_ 

Y 


Thus , 

1 12 
Soo  = y ( ^ ~ 1 

The  turbulent  viscosity  coefficient  and  the  source  term  are 
set  to  TURMU(L,J)  = 0 and  S(L,N,J)  = 0,  respectively,  where  N 
= 1~5  represents  five  fluid  properties.  There  is  an  "IF" 
statement  controlled  by  IREAD  such  that 


IREAD 


0 ; no  reading  and  passing  over  this  statement 
for  calculating  the  initial  conditions 


1 ; reading  a new  set  of  input  for  restart 


The  last  process  is  to  scale  the  variables  by  Jacobian 


A A A A 


(4.26) 


where  £ = 1,  . . .,  LMAX ; n = 1,  . . .,  5;  and  j = 1 , . . ., 
JMAX. 

SUBROUTINE  EIGEN 

This  subroutine  calculates  the  maximum  Courant  number, 
with  CNBR  = 0.0  as  input,  and  the  At  with  CNBR  ^ 0.0  as 
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input.  The  purpose  of  this  subroutine  is  to  provide  the 
information  for  the  stability  of  the  calculation. 

SUBROUTINE  STEP 

An  implicit  integration  scheme,  Eqn.  (4.24),  is 
implemented  by  this  subprogram.  There  are  six  subroutines, 
BC,  RHS , SMOOTH,  FILTRX,  FILTRZ , and  BTRI , called  here.  The 
step  algorithm  performs  the  computations  described  in  order 
as  follows: 

(1)  CALL  SUBROUTINE  BC 

The  subprogram  BC  calculates  the  fluid  properties  on  the 
boundary.  Since  the  outer  boundary  is  in  the  free  stream 
region,  the  values  should  remain  in  the  free  stream 
condition.  For  the  inner  boundary  condition,  the  velocities 
are  calculated  by  Eqn.  (4.13),  the  body  surface  pressures 
obtained  by  the  implementation  of  Eqn.  (4.18),  and  the  total 
flow  energies  per  unit  volume  along  body  surface  are  found  by 
Eqn.  (3.8).  The  fluid  properties  on  the  axis  boundary  (AB 
line  in  Fig.  2.1)  are  obtained  by  making  an  extrapolation 
from  interior  points  of  the  domain,  and  the  same  is  done 
along  the  downstream  boundary  (DE  line  in  Fig.  2.1).  The 
values  on  the  four  segments  of  the  boundary  are,  therefore, 
determined. 

(2)  CALL  SUBROUTINE  RHS 

This  subprogram  constructs  and  calculates  the  vector 
form  S']  of  source  term,  one  part  of  the  right  hand  side  of 


Eqn.  (4.22), 
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51  = - At[ 6£En+6?Gn-Re-1 5^Sn]  - AtHn 

(3)  CALL  SUBROUTINE  SMOOTH 

In  the  other  part  of  the  right  hand  side  of  Eqn.  (4.22), 
the  explicit  smoothing  terms  S£ 

52  = - eEJ-1 [(A5V5)2+(A^V?)2]Jqn 
are  executed. 

(4)  Once  we  have  obtained  the  values  of  the  source 
terms  in  Eqn.  (4.22),  the  residual  at  the  end  of  10  time-step 
is  computed  by  taking  the  root  mean  square  of  these  summation 
values,  such  that 


S 


JMAX, 5 , LMAX 

I 

j=k=£=1 


Residual 


S 

( JMAX-2 ) ( LMAX-2 ) 


and 


Thus,  the  residual  is  the  mean  residual  over  the  whole 
field.  This  value  is  printed  out  for  the  inspection  that 
gives  some  information  on  the  accuracy  of  the  results  for  the 
purpose  of  use. 

(5)  CALL  SUBROUTINE  FILTRX 

This  subprogram  constructs  and  calculates  the 
coefficients  of  the  unknown  vector  Aq*  j.n  Eqn.  (4.24(a));  the 
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set  of  these  coefficients  is  formed  by  a block  tri-diagonal 
matrix  with  5-factor, 

[^-matrix]11  = [I+h5  An-eTJ-1  V A J]  , l = 2,  3,  . . . , LM 

(6)  CALL  SUBROUTINE  BTRI 

Thus,  Aq*  is  found  calling  this  subprogram  which  is  the 
block  tri-diagonal  solver. 

(7)  CALL  SUBROUTINE  FILTRZ 

A ^ 

Since  Aq*  becomes  known  vector,  we  then  substitute  Aq* 
into  Eqn.  (4.24(b)).  Similarly,  this  subprogram  is  called  to 
evaluate  a block  tri-diagonal  5-matrix,  the  coefficients  of 
the  unknown  vector  Aq  in  Eqn.  (4.24(b))-,  formed  as 

[5-matrix]n  = [I  + hS  Cn-eTJ  ^A_V_J-hRe  ^ <5  J "^Mnj]  , 

5 I C 5 ? 

j = 2,  . . . , JM 

(8)  CALL  SUBROUTINE  BTRI 

The  unknown  vector  Aqn  is  found  by  calling  BTRI  once 
more.  Thus,  the  qn+1  in  Eqn.  (4.24(c))  can  be  readily  found. 
SUBROUTINE  MAP 

This  subroutine  sketches  the  digital  map  of  the 
calculated  pressure  distribution  for  assistance  to  inspect 
the  behavior  of  pressure  calculation  roughly. 
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SUBROUTINE  OUTPUT 

The  subroutine  OUTPUT  is  designed  to  print  out  the  dis- 
tribution of  fluid  properties  in  the  domain  or  on  the 
boundary. 


CHAPTER  V 

ADAPTIVE  GRID  GENERATION  CODE 

The  geometry  of  a secant-ogive-cylinder-boattail  (SOCBT) 
projectile  is  shown  in  Fig.  5.1;  the  model  has  a 5-caliber 
secant-ogive,  a 2-caliber  cylinder,  and  a 1-caliber  7° 
boattail.  Since  the  projectile  is  axisymmetric , we  just  need 
one-half  the  entire  domain  as  shown  in  Fig.  2.1.  The  outer 
boundary  is  selected  to  keep  free  stream  everywhere  on  this 
boundary.  Behind  the  boattail  there  exists  the  wake  flow 
which  is  not  easy  to  solve.  It  is  modeled  as  an  attached 
"sting"  instead  of  the  wake  for  simplicity.  The  axis 
boundary  (5  = 1,  AB  line  in  Fig.  2.1)  coincides  with  the 
model  axis  and  the  downstream  boundary  (5  = M,DE  line  in  Fig. 
2.1)  is  perpendicular  to  the  sting. 

To  illustrate  the  use  of  the  code  GRIDGEN  mentioned  in 
Chapter  I,  it  is  applied  to  the  projectile  flow  problem  for 
numerical  experiment  in  this  study.  For  example,  the  flow 
field  has  the  computational  domain  extended  to  4-times  body 
length  in  front,  5 body  lengths  between  the  projectile  axis 
and  the  straight  line  of  outer  boundary,  and  5 body  lengths 
behind  the  actual  projectile.  The  sting  is  modeled  (see  Fig. 
2.1)  as  a natural  extension  of  the  boattail  for  a distance  of 
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Fig.  5.1  Projectile  configuration 
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1.767  calibers  and  then  turned  parallel  to  the  model  axis. 
Thus,  the  boundaries  of  the  flow  field  are  formed.  There  are 
70x35  grid  points  used  in  this  domain,  70  points  along  the  5- 
direction  (23  points  on  the  ogive,  22  points  on  the  cylinder, 
17  points  on  the  boattail,  and  8 points  on  the  sting),  and  35 
points  along  the  c-direction.  The  grid  points  of  the  inner 
and  outer  boundaries  are  defined  initially  by  using  the 
geometrical  analysis  techniques  for  the  projectile  flow 
computation  developed  by  Steger,  Nietubicz,  and  Heavey 
[16].  Then,  the  grid  points  on  the  boundary  are  regenerated 
adaptively  with  a cubic  function's  clustering  routine  which 
permits  the  user  to  choose  regions  of  clustering  according  to 
the  user's  experience  and  intuition.  The  grids  in  the  first 
subdomain,  the  ogive  part,  are  generated  by  using  an  elliptic 
solver;  the  initial  locations  of  grid  points  in  this 
subdomain  are  specified  by  equal  spacing  interpolation  on  the 
lines  which  are  the  connections  between  the  inner  boundary 
point  and  the  corresponding  outer  boundary  points  [see  the 
dot  points  in  the  first  subdomain  of  Fig.  5.2(a)].  The  other 
parts,  cylinder,  boattail,  and  sting,  are  the  second 
subdomain  which  can  be  readily  generated  by  the  straight  ray 
technique  with  the  clustered  grid  points  on  the  boundary  as 
shown  in  the  second  subdomain  in  Fig.  5.2(a).  The  grid 
network  is  then  generated;  also,  this  grid  network  can  be 
modified  to  emphasize  high  grid  resolution  in  the  viscous 
layer  by  clustering  grid  points  with  an  exponential  function 
along  grid  lines  normal  to  the  streamwise  direction  as  shown 
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(a)  Grid  points  clustering  on  the  inner  and  outer  boundaries 
by  GRIDGEN 


(b)  Grid  points  clustering  close  to  the  inner  boundary  by 
GRIDGEN 


Fig.  5.2  Grid  networks  generated  by  GRIDGEN 
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(c)  Grid  network  generated  by  GRIDGEN 
Fig.  5 .2-continued. 


Fig.  5-3  Initial  grid  network  generated  by  adaptive  GRIDGEN 


87 


in  Fig.  5.2(b).  The  resultant  grid  network  generated  by 
GRIDGEN  is  shown  in  Fig.  5.2(c). 

As  mentioned  in  the  Introduction,  the  current  version  of 
GRIDGEN  can  generate  a planar  grid  network  by  solving  either 
an  elliptic  boundary  value  problem  or  a hyperbolic  initial 
value  problem.  Also,  Hsu  [29]  has  provided  an  additional 
option  to  GRIDGEN  for  generating  a grid  network  with  the 
resolution  of  boundary  grid  adaptive  to  a specified  control 
function  such  as  a computed  pressure  gradient  distribution 
instead  of  boundary  grid  adaptive  to  user's  intuition. 
Moreover,  he  modified  the  hyperbolic  grid  routines  of  GRIDGEN 
for  better  smoothness  and  orthogonality  of  grid  networks. 

Using  the  same  size  of  domain  and  the  same  number  of 
corresponding  grid  points  on  the  boundaries  as  used  in 
GRIDGEN,  the  initial  grid  network  adaptive  to  the  control 
weight  function  is  generated  by  using  N-R  method  with  SOR(a)  = 
1.6)  to  solve  Eqn.  (2.57)  within  the  entire  domain.  Since  no 
information  is  available  initially  for  measuring  the 
concentration  of  the  grid,  we  set  X = x/x^  = 1,  Xw  = X / X w = 
0.  The  grid  points  on  the  inner  and  outer  boundary  are 
specified  as  either  the  uniform  or  the  arbitrary  grid 
spacing.  Since  experience  shows  that  the  grids  on  the  axis 
(AB  line  in  Fig.  5.3)  and  downstream  (DE  line)  boundaries 
adaptive  to  the  control  function  obtained  from  the  solution 
of  the  thin-layer  code  do  not  give  a good  convergent 
solution,  these  grid  points  are  determined  by  an 
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extrapolation  of  the  interior  grid  curves  to  be  orthogonal  to 
the  boundary.  This  extrapolating  procedure  is  similar  to 
that  used  to  obtain  the  fluid  properties  on  the  boundary  with 
the  thin-layer  code  in  Chapter  IV.  The  iteration  continues 
until  either  the  solutions  converge  everywhere  to  within  the 
allowable  error  or  it  is  quitted  at  a specified  maximum 
iteration  number.  The  resulting  initial  grid  network,  for 
example,  is  shown  in  Fig.  5.3. 

Once  the  approximate  solution  for  the  fluid  properties 
is  obtained  from  the  calculation  of  50-100  iterations  by  the 
thin-layer  code  based  on  the  initial  grid  network;  the  grid 
generations  adaptive  to  some  property  of  the  given  solution, 
which  is  the  pressure  gradient  used  in  this  code  as  the 
weight  function,  are  implemented  with  the  scaled  Aw  and  XQ  to 
solve  Eqn.  (2.64)  and  (2.57)  on  the  boundary  and  within  the 
domain,  respectively. 

For  successful  adaptive  grid  network  generation,  a 
weight  function  is  needed  that  is  a measure  of  the  steepness 
of  the  gradient  of  the  important  physical  phenomenon.  For 
example,  a steep  pressure  gradient  may  occur  at  the 
corresponding  location  where  a shock  wave  front  exists  in  the 
flow  field  of  the  high  speed  aerodynamics  problem.  Thus,  a 
norm  of  the  pressure  gradient,  equal  to  the  sum  of  the 
absolute  values  of  the  Cartesian  components  of  the  pressure 
gradient  is  used  in  the  solutions  given  here,  such  that 
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W(x,z)  = | Px | + |Pz|  in  n 


(5.1) 


V(s)  = |Ps| 


on  r 


(5.2) 


where  Px  = P?CX+P??X*  pz  = p^z+p^z  are  found  by  applying 


Eqn.  (2.27)  to  the  computational  form 


(5.3) 


The  Ps  in  Eqn.  (5.2)  is  the  absolute  value  of  pressure 
gradient  along  the  boundary,  such  that  Ps  = P^/S^. 

Since  we  discretize  W(x,z)  in  central  differences,  it 
must  be  noted  that  the  pressure  gradient  at  peak  points 
(i.e.,  maximum  and  minimum  value  of  pressure)  causes  a large 
truncation  error.  Thus,  the  adaptive  grid  may  not  generate 
an  optimum  grid  point  in  this  region;  the  pressure  gradient 
at  certain  peak  points  may  be  treated  specially  by  taking 
averages  in  the  vicinity  of  this  point,  such  that 


Also,  because  the  finite  difference  equation  roughens  the 
data  causing  spurious  changes  for  W(x,z),  Brackbill  and 
Saltzman  suggested  the  smoothing  equations, 


(5.4) 
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Wj,  k 


(5.5) 


W£+2 

o»k 


= w* 


j»k 


+ v[(W 


+w^ 

j+1,k+wj-1,k 


)/2_WS,k]  on  r 


(5.6) 


"to  modify  the  weight  functions.  These  equations  are  executed 
with  v = 0.4  for  two  or  three  iterations,  which  produces  a 
good  result. 

To  solve  the  Eqn.  (2.48),  the  boundary  grid  system  must 
be  transformed  from  2-D  to  1-D.  The  line  space,  S = S(x,z), 
is  defined  as 


S1  = 


M 


1 

2 n 2 


h [(xj+1,1"xj,l)  +(zj+1,l"zj,l)  ] 


along  inner  r 


(5.7) 


M 


sd  ■ 


1 

2,2 


[(X0+1 ,N"xj,N}  +(zj+1 ,N"zj,N}  ] 


along  outer  r 


(5.8) 


For  inner  r,  there  are  four  segments  (ogive,  cylinder, 
boattail,  and  sting)  combined  together,  the  superscript  i is 
the  number  of  the  segment,  for  example,  Sj  is  the  line  space 
on  the  cylinder.  This  treatment  is  required  to  fix  the  two 
end  points  of  each  segment  in  order  to  keep  the  projectile 
shape  unvarying  when  implementing  the  iterations  for  boundary 
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adaptive  grid  generation.  In  other  words,  the  adaptive  grid 
for  the  inner  boundary  is  executed  separately,  segment  by 
segment.  The  S-j  on  outer  r is  implemented  without  any 
separation.  Once  we  have  solved  for  we  go  back  to  find 
the  corresponding  (x^-pZj^)  on  the  inner  boundary  and 
(x  j , N> z j , Ij)  on  the  outer  boundary  by  employing  the  cubic 
spline  interpolation  function. 

The  adaptive  grid  generation  code  with  control  weight 
function  is  primarily  based  on  the  code  GRIDGEN.  It  is 
modified  by  the  insertion  of  a few  subroutines  and  provides  a 
new  option.  The  details  of  the  new  subroutines  will  be  found 
in  Appendix  B.  Also,  it  is  convenient  by  coupling  the 
adaptive  GRIDGEN  with  control  function  and  thin-layer  code 
together  to  obtain  a direct  computation  for  projectile  flow 
problem.  One  can  then  change  the  input,  the  new  grid,  and 
the  time  level  automatically  with  a few  control  numbers. 


CHAPTER  VI 

NUMERICAL  EXPERIMENTS 


As  mentioned  in  Chapter  V,  the  thin-layer  code  and 
adaptive  grid  generation  code  are  coupled  together  for 
calculating  the  solution  of  the  projectile  flow  problem.  The 
thin-layer  code  provides  the  prior  approximate  solution 
properties  for  the  control  weight  function  of  the  adaptive 
grid  generation  code.  This  grid  generation  code  then 
provides  the  adaptive  grid  network  for  the  thin-layer  code 
calculation  of  fluid  properties.  This  cycle  is  repeated 
until  a satisfactory  solution  is  found. 

For  assessing  the  effectiveness  of  the  adaptive  grid 
generation  technique  developed,  we  have  considered  the 
axisymmetric  inviscid  transonic  projectile  aerodynamics 
problems.  The  smoothing  parameters  are  chosen  as 

= 2£g  = 8At  for  this  study.  The  weight  function  employed 
for  grid  resolution  can  be  related  to  the  solution  of  fluid 
properties  such  as  the  pressure,  velocities,  density  and 
temperature.  However,  the  pressure  gradient  seems  to  be,  in 
some  sense,  a rather  strong  gradient  of  physical  phenomena 
describing  the  location  of  shocks  in  the  inviscid  fluid  flow. 

Brackbill  and  Saltzman  chose  the  high  power  of  the 
physical  gradient  as  the  weight  function,  such  as  in 
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reference  [8]  for  solving  inviscid  supersonic  problems  and 
(—rfr) ^ in  reference  [7]  for  a steady  convection  diffusion 
problem.  For  the  choice  of  a proper  weight  function  in  the 
projectile  flow  field,  we  feel  that  the  first  order  power  of 
Vp  will  be  sufficient  enough  for  emphasizing  the  concentra- 
tion of  grid.  Therefore,  the  weight  function  used  in  this 
study  is  defined  as 

W(x,z)  = |PJ  + |P2|  . 

For  conducting  the  numerical  investigations  on  the 
adaptive  grid  generation  technique,  an  experimental  case  for 
the  calculation  with  the  weight  function,  W(x,z)  = |PX|  + 
|PZ|,  and  1M  = 0.91  was  executed  by  the  coupled  computer 
codes.  The  constant  penalty  parameters  *w  = and  x0  =' 

VXq  defined  in  Eqn.  (2.20)  were  used  in  this  case,  in  which 
and  are  constants  similar  to  those  of  Saltzman  which 
have  been  described  in  Chapter  II,  and  * = 4.0.  The  initial 
grid  network  used  is  the  one  shown  in  Fig.  5.3  (70x35 
points).  The  computed  surface  pressure  coefficient  Cp  at  NT 
= 50  (NT  stands  for  the  number  of  time  steps)  is  shown  in 
Fig.  6.1(a);  also,  the  measured  data  (diamonds)  for  steady 
flow  obtained  from  reference  [30]  are  plotted  for  comparison 
with  the  solution.  A new  grid  network  generated  adaptively 
to  the  approximate  solution  at  NT  = 50  is  shown  in  Fig. 
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X/D 


Fig.  6.1  Surface  pressure  and  adaptive  grid  networks' 
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6.1(b).  Subsequently,  a new  updated  adaptive  grid  network  is 
generated  at  every  150  time-step  intervals.  The  computa- 
tional process  proceeds  to  NT  = 350;  however,  in  the  next  150 
time  steps,  the  thin-layer  Navier-Stokes  code  fails  to  con- 
verge. Figure  6.1(a)  shows  the  surface  pressure  distribution 
for  NT  = 50,  which  tends  to  converge  to  the  measured  data. 

The  grid  network  generated  adaptively  to  the  approximate 
solution  at  NT  = 50  gives  a large  grid  spacing  near  the 
region  of  the  boattail  as  shown  in  Fig.  6.1(b).  The  next  150 
time-step  calculation  gives  the  Cp-distribution  at  NT  = 200 
in  Fig.  6.1(c).  Apparently,  the  pressure  near  the  large  grid 
spacing  region  in  Fig.  6.1(b)  does  not  agree  with  the 
measured  data.  In  fact,  it  has  a positive  peak  where  the 
measured  data  have  a negative  peak.  These  poor  results 
probably  come  from  the  large  grid  spacing  described  in  Fig. 
6.1(b).  Similarly,  the  grid  spacing  generated  adaptively  to 
the  solution  at  NT  = 200  is  also  large  near  the  boattail 
region  in  Fig.  6.2(a).  The  surface  pressure  calculated  at  NT 
= 350  disagrees  greatly  with  the  measured  data  near  that 
region  in  Fig.  6.2(b).  Consequently,  the  grid  spacing  of  the 
grid  network  generated  adaptively  to  the  solution  at  NT  = 350 
is  still  large  in  the  same  region  as  shown  in  Fig.  6.2(c). 

At  last,  these  new  grid  networks  generated  adaptively  to  an 
unreasonable  solution  cause  the  computational  instability. 
This  large  grid  spacing  may  be  caused  by  the  smoothing 
effect,  as  occurred  in  conventional  elliptic  grids,  and  the 
small  value  of  the  control  weight  function  in  that  region. 
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(a)  Grid  network  generated  at  NT  = 200 


Surface  pressure  at  NT  = 350 


(c)  Grid  network  generated  at  NT  = 350 


Fig.  6.2  Surface  pressure  and  adaptive  grid  networks' 
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To  overcome  this  difficulty,  an  exponential  clustering 
similar  to  that  of  GRIDGEN  is  applied  to  all  grid  lines 
normal  to  the  streamwise  direction.  The  first  grid  spacing 
at  the  boundary  is  specified  as  aS  = 0.01.  The  resulting 
adaptive  grid  network  with  clustering  is  shown  in  Fig.  6.3(a) 
and  the  effect  of  the  clustering  can  be  seen  by  comparing 
Fig.  6.3(b)  and  Fig.  6.2(c).  The  application  of  this 
clustered  grid  network  and  the  use  of  the  poor  results 
computed  at  NT  = 350  to  restart  the  computational  process 
results  in  stable  computation.  The  surface  pressure 
coefficient  computed  after  150  time  steps  is  clearly 
converging  to  the  measured  data  as  shown  in  Fig.  6.3(c).  To 
confirm  the  importance  of  a clustering  strategy,  a new  grid 
is  generated  at  NT  = 500  for  the  continuation  of  the 
integration  process.  The  results  obtained  at  NT  = 650  and 
shown  in  Fig.  6.4(b)  indicate  that  the  implementation  of  the 
clustering  strategy  in  the  adaptive  grid  generation  can 
provide  a good  grid  network  for  the  projectile  aerodynamics 
computations . 

The  adaptive  grid  generation  with  clustering  was  next 
tested  on  the  projectile  aerodynamic  problem  at  1M  = 0.96  to 
see  whether  a better  grid  network  was  obtained,  which  would 
give  better  results.  The  adaptive  grid  system  is  changed 
after  each  200  time-step  interval.  The  constant  penalty 
parameters  and  weight  function  used  are  the  same  as  in  the 
previous  experimental  case.  The  Cp-distr ibution  computed  at 
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(a)  Grid  network  with  clustering  ( AS  = 0.01)  at  NT  = 350 


(b)  Expanded  view  of  grids  near  the  body  surface 


(c)  Surface  pressure  calculated  on  the  above  grids  at  NT  = 500 

Fig.  6.3  Adaptive  grid  networks  and  surface  pressure 
distribution 


99 


(a)  Grid  network  generated  at  NT  = 500  with  clustering 


(b)  Surface  pressure  at  NT  = 650 


Fig.  6.4  Adaptive  grid  network  and  surface  pressure 
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NT  = 1650  is  plotted  in  Fig.  6.5(a).  It  is  clear  that  some 
values  of  the  surface  pressure  Cp  still  do  not  approach 
closely  to  their  corresponding  experimental  data. 

Apparently,  we  still  cannot  obtain  a better  adaptive  grid 
network.  In  order  to  improve  the  adaptive  grid  network 
further,  Brackbill  and  Saltzman  paid  much  attention  to 
modifying  their  weight  function  and  to  getting  another  type 
of  weight  function  [31].  Also,  high  order  power  of  physical 
gradient  was  chosen  as  weight  function  in  their  reports,  as 
described  previously.  All  of  those  treatments  for  the 
weight  function  were  in  order  to  emphasize  the  concentration 
of  grids  more  strongly. 

The  most  fruitful  direction  for  obtaining  a good  grid 
system  is  probably  in  the  determination  of  the  penalty 
parameters.  Unfortunately,  the  value  of  X ^ = W(L/L')^ 
defined  in  Eqn.  (2.24)  within  the  domain  is  difficult  to 
match  to  the  values  of  X^  = W(L/L')^  given  in  Eqn.  (2.45)  on 
the  boundaries  when  we  apply  these  normalized  constant 
penalty  parameters  to  the  1-D  variational  principle  for  the 
adaptive  boundary  grid.  Since  the  four  side  segments  of  the 
boundary  have  different  lengths  of  L corresponding  to  the 
number  of  grid  points  L' , the  values  of  (L/L1)  for  these 
segments  are  not  equal.  Consequently,  to  gain  the 
consistency  between  boundary  and  domain,  the  value  of  (L/L') 
in  the  domain  cannot  be  constant  in  the  domain.  Otherwise, 
if  a reasonable  value  of  (L/L')>  such  as  the  square  root  of 
the  total  grid  area  over  the  total  grid  points,  is  chosen 
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(a)  Body  surface  pressure  at  NT  = 1650  with  constant 
penalty  numbers 


(b)  The  variance  over  the  mean  as  a function  of  Xw 
Fig.  6.5  Surface  pressure  and  the  variance  versus  Xw 
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for  the  domain  the  investigation  of  an  optimal  value  of  the 
global  parameter  X in  Xw  = X/X^  becomes  tedious,  and  it  may 
not  give  a good  result.  Therefore,  the  constant  penalty 
parameters  may  not  be  suitable  for  generating  a good 
adaptive  grid  network  in  both  boundary  and  domain.  The  con- 
sistency of  the  grid  points  between  boundary  and  domain  thus 
becomes  extremely  important  to  obtain  the  accurate  results. 

In  order  to  demonstrate  substantially  the  properties  of 
WJ  = const,  in  the  projectile  flow  field,  an  experiment 
was  performed  to  investigate  the  behavior  of  the  variation 
of  the  concentration  functional  versus  X^.  jn  this  experi- 
ment, XQ  is  set  to  1,  and  Xw  varies.  The  pressure  gradient 
calculated  at  NT  = 1650  is  chosen  as  the  weight  function. 

We  then  use  statistical  analysis  to  find  the  dependence  on 
Xw  of  some  measures  of  the  distribution  of  WJ-2.  The  mean 
and  variance  of  the  variation  of  WJ-2  are  defined  as 


mean 


M,N 

[ l 

j=k=1 


(WJ-2).  , ]/MN 

J 9 & 


M,N 

variance  = [ l (WJ  -mean)  . ]/MN 
j=k=1 


(6.1) 

(6.2) 


Table  6.1  shows  the  variance/mean  of  the  values  WJ~2  as  X 

w 

varies.  These  data  are  also  plotted  in  Fig.  6.5(b)  which 
clearly  indicate  that  the  variance  decreases  and  approaches 
to  the  solution  WJ-2  = const,  as  xw  increases.  These 
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Table  6.1 

Numerical  Experiment  Data  for  Variance/Mean  Versus  Xw 


x 


w 


Variance/Mean 


0.00001 
0.00005 
0.00010 
0.00050 
0.00100 
0.00500 
0.01000 
0.05000 
0.10000 
0.50000 
1 .00000 
2.00000 
4.00000 
8.00000 
10.00000 


0.9034E+01 
0.9034E+01 
0.9034E+01 
0.9032E+01 
0.9030E+01 
0.9013E+01 
0.8993E+01 
0.8845E+01 
0.8704E+01 
0.8279E+01 
0.81 39E+01 
0.8032E+01 
0. 7943E+01 
0.7863E+01 
0.7839E+01 


results  may  be  helpful  for  finding  adequate  penalty 
parameters  later  on. 

As  discussed  previously,  X^  and  X^  may  be  considered  as 
the  variables  within  the  domain  and  on  the  boundaries  in 
order  to  emphasize  the  orthogonality  and  concentration 
locally  and  obtain  a good  consistency  of  grids  between 
boundary  and  domain.  Thus,  the  calculations  of  L,  L'  and  W 
may  be  redefined  in  each  cell,  such  that 
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L = aS(  5 » 0 = < 


|D||  = [(X  )2+(Z  )2]2 

j,k  5 5 3,k 


|Dt|  = [(X  )2+(Z  )2]2 

5 5 j,k 


along  5-  (6.3) 

direction 


along  5-  (6.4) 

direction 


L*  = 1,  and  ¥ = W(5>c)  (6.5) 

where  l^k  k and  |Dc  | . , are  the  magnitudes  of  the  vectors, 
{D?}  and  {D^},  respectively,  defined  in  Eqn.  (2.12)  at  point 
(j,k)  in  the  domain  and  on  the  boundary.  Thus,  L becomes 
the  grid  spacing  aSj ^ and  W represents  the  corresponding 
value  of  weight  function  in  the  computational  space. 

Consequently,  the  parameters  x^  and  x£  in  xw  = x/x^  and  x0  = 
x/Xq,  respectively,  are  converted  into  variables  given  by 


(x0}j,k  - 

Xo(5>c) 

^xw^j,k 

(xw} j,1  " 

5 > C ) 

^ Xw  ^ j , N = 

*iU»c) 

(is4)d,k 


Wd>k(4S4) 
w,  ^AS3) 

J 9 1 


¥.>n(aS") 


, in 

d,1  ’ 

j,N  ’ 


ft 

in  a 

on  inner 
on  outer 


r 


r 


► (6.6) 


(6.7) 


The  global  parameter  x may  be  specified  to  depend  on  At,  and 
will  be  assumed  as  x = lOOAt  in  an  example.  After 
substitution  of  Eqn.  (6.7)  and  x into  Eqn.  (2.59),  the  grid 
resolutions  on  the  boundary  are  readily  obtained  by 
executing  the  Eqn.  (2.64).  Once  the  boundary  conditions  are 
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known,  we  then  substitute  Eqn.  (6.6)  and  X into  (2.51)  to 
obtain  the  grid  resolutions  within  the  domain  by  Eqn.  (2.57) 
in  alternating  directions,  such  that 


Iterations  = Odd  numbers  Iterations  = Even  numbers 


J=2, 3 , . . . ,M-1 

i 1 

| Eqn.  (2.57)  with  [ 

Eqn . (6.4) 

! K=2,3, • • • ,N-1 

I I 


K=2 , 3 , . . . , N-1 

I 

Eqn.  (2.57)  with 
Eqn.  (6.3) 

| J=2, 3 , • • • >M-1 


The  adaptive  grid  generation  based  on  the  variation 
principles  but  implemented  with  a clustering  strategy  and 
variable  penalty  parameters  X^  and  was  investigated  on 
the  projectile  flow  problem  at  1M  = 0.91,  0.96  and  1.10. 

The  local  penalty  parameters  and  were  tested  in  those 
investigations  with  the  expectation  of  obtaining  good 
results.  The  strategy  of  this  entire  computation  for  the 
projectile  inviscid  flow  problem  is  shown  in  Table  6.2, 
where  NG  is  the  number  of  changes  of  new  grid  system,  DT  is 
the  time  step  At,  X is  the  global  penalty  parameter,  NT  is 
the  number  of  time  steps  and  DS  is  the  clustering  first  grid 
spacing  right  next  to  the  boundary  surface.  A sequence  of 
the  computed  surface  pressure  coefficients  Cp  and  the 
generated  adaptive  grid  network  for  the  case  of  free  stream 
Mach  number  1M  = 1.10  is  shown  in  Figs.  6. 6-6. 9.  This 
sequence  indicates  clearly  the  process  by  which  the  grid 
network  generated  is  adaptive  extremely  well  to  the  pressure 
gradient  field  and  exhibit  the  shock  pattern.  The  computed 
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Table  6.2 

Strategy  of  Implementation  for  Each  NT  = 150  Interval 


NG 

X 

DT 

NT 

DS 

0 

GENERATE  THE 

ORIGINAL  GRID 

0.01 

COMPUTE 

0.005 

1-50 

1 

0.5 

CHANGE  TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.005 

0.010 

0.010 

51-100 

101-150 

151-200 

2 

1 .0 

CHANGE  TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.025 

0.025 

0.050 

201-250 

251-300 

301-350 

3 

5.0 

CHANGE  TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

0.075 

0.075 

351-400 

401-450 

451-500 

4 

7.5 

CHANGE  TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.100 

0.100 

0.100 

501-550 

551-600 

601-650 

5 

10.0 

CHANGE  TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.100 

0.100 

0.100 

651-700 

701-750 

751-300 
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GRID  NETWORK  GENERATED  AT  NT=  50 


Fig.  6.6  Surface  pressure  and  resulting  grid  network  at 
NT  = 50 
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GRID  tETUORK  GENERATED  AT  NT=  200 


Fig.  6.7  Surface  pressure  and  resulting  grid  network  at 
NT  = 200 
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GRID  NETWORK  GENERATED  AT  NT=  350 


Fig.  6.8  Surface  pressure  and  resulting  grid  network  at 
NT  = 350 
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GRID  NETWORK  GENERATED  AT  NT=  5M 


Fig.  6.9  Surface  pressure  and  resulting  grid  network  at 
NT  = 500 
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surface  pressure  coefficient  Cp,  for  1M  = 0.91 , 0.96,  and 
1.10  obtained  at  NT  = 650,  800,  and  500,  respectively,  are 
shown  in  Fig.  6.10.  Each  of  the  three  flow  cases  has  been 
solved  again  to  NT  = 800  using  a fixed  but  good  adaptive 
grid  system  generated  from  the  original  GRIDGEN  as  described 
in  Chapter  V and  the  results  obtained  are  shown  in  Fig.  6.11 
for  comparison.  It  is  clear  that  the  results  obtained  from 
adaptive  gridding  agree  better  with  the  experimental  data. 

Additional  numerical  experiments  have  been  designed  and 
conducted  to  gain  some  insight  into  the  relative  importance 
of  adaptive  boundary  grids  and  internal  grids  to  the 
accuracy  and  solution  algorithm  of  the  Navier-Stokes  code. 
For  each  of  the  flow  cases  considered,  a planar  grid  network 
is  generated  from  the  original  GRIDGEN  with  the  surface 
boundary  grid  distribution  exactly  the  same  as  that  of  the 
corresponding  adaptive  grid  network.  Figure  6.12  shows  the 
difference  of  the  two  grid  systems  for  1M  = 1.10.  The  grid 
network  generated  is  then  employed  in  the  Navier-Stokes  code 
for  the  computation.  Accurate  results  have  been  obtained 
for  the  case  1M  = 1.10;  however,  for  the  other  two  cases  1M 
= 0.91  and  0.96,  the  solution  algorithm  of  the  Navier-Stokes 
code  failed  within  NT  = 30.  Thus,  for  the  fixed  grid 
network,  the  original  GRIDGEN  may  not  generate  a good 
interior  grid  network  within  the  domain  based  on  an  adaptive 
boundary  grid. 

In  order  to  find  a better  adaptive  grid  network 
further,  there  may  be  a better  strategy  implemented  by 
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Fig.  6.10  Converged  solutions  (CD)  obtained  from  adaptive 
grid  system 
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Fig.  6.11  Body  surface  pressure  calculated  on  fixed  grid 
network  at  NT  = 800 
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(a)  Adaptive  grid  network  generated  at  NT  = 500 


(b)  An  adaptive  grid  system  obtained  from  the  original 

GRIDGEN  with  exactly  the  same  surface  boundary  grids  as 
those  on  Fig.  6.12(a) 

Fig.  6.12  Adaptive  grid  networks  generated  for  1M  = 1.10 
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reducing  the  final  time  step  from  At  = 0.1  to  0.5  and 
increasing  Un  ^ = X/Aw  and  Xo  = X^X'o  gradually  as 
described  in  Table  6.3.  Figures  6.13  and  6.14  show  the 
adaptive  grid  networks  generated  for  each  100  time  steps  in 
the  case  1M  = 0.96  and  the  converged  Cp-distribution  calcu- 
lated at  NT  = 1000.  These  results  agree  more  smoothly  and 
better  with  the  experimental  data  than  the  results  shown  in 
Fig.  6.10  for  1M  = 0.96.  We  then  generate  a planar  grid 
network  by  the  original  GRIDGEN  with  the  surface  boundary 
grid  distribution  exactly  the  same  as  the  corresponding 
adaptive  grid  network  in  Fig.  6.13;  this  fixed  grid  system 
is  employed  in  the  thin-layer  Navier-Stokes  code  for  com- 
putation and  provides  the  results  as  shown  in  Fig.  6.15.  It 
is  clear  that  the  results  obtained  from  the  adaptive  grid  in 
Fig.  6.13  agree  better  with  the  experimental  data  than  the 
results  obtained  from  the  fixed  grid  in  Fig.  6.15. 

Also,  from  the  results  shown  in  Fig.  6.15,  the 
calculated  surface  pressure  Cp  on  the  ogive  and  cylinder 
agrees  well  with  the  experimental  data;  however,  the  surface 
pressure  is  far  away  from  the  measured  data  between  the 
boattail  and  sting.  This  is  caused  by  a grid  network  being 
generated  which  is  not  adaptive  to  the  solution.  Thus,  we 
conclude  that  satisfactory  results  depend  strongly  not  only 
upon  the  adaptive  grid  points  on  the  boundary  but  also  upon 
the  adaptive  grids  within  the  domain. 

For  summarizing  the  strategies  of  Tables  6.2  and  6.3, 
we  also  conclude  that  the  global  penalty  parameters  * cannot 
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Table  6.3 

Strategy  of  Implementation  for  Each  NT  = 100  Interval 


NG 

X 

DT 

NT 

DS 

0 

GENERATE  THE  ORIGINAL  GRID 

0.01 

COMPUTE 

0.005 

1-100 

1 

1 .0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.010 

101-200 

2 

1.5 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.025 

201-300 

3 

3.0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

301-400 

4 

6 . 0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

401-500 

5 

10.0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

501-600 

6 

16.0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

601-700 

7 

23.0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

701-800 

8 

28.0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

801-900 

9 

30.0 

CHANGE 

TO 

THE 

NEW  GRID 

0.01 

COMPUTE 

0.050 

901-1000 
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GRID  tCTUORK  GENERATED  AT  NT=  9M 


Fig.  6.13  Adaptive  grid  network  and  converged  solution  (Cp) 
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(a)  Adaptive  grid  network  generated  at  NT  = 1000 


(b)  The  grid  s-line  exhibiting  the  shock  patterns 

Fig.  6.14  Grid  network  adaptive  to  the  converged  solution 
at  NT  = 1000 
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FIXED  GRID  NETWORK  GENERATED  BY  "CRIDGEN* 


Fig.  6.15  Adaptive  grid  system  obtained  from  "GRIDGEN"  with 
exactly  the  same  surface  boundary  grids  as  those 
on  Fig.  6.15  and  its  results  calculated  at  NT  = 
1000 
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be  used  with  a large  number  initially  in  the  process  of 
calculation;  however,  X can  be  increased  when  the  final 
time-step  is  fixed  in  the  continuing  computation. 

Therefore,  the  penalty  parameters  xw  and  XQ  can  be  certainly 
converted  into  variables  which  are  not  only  in  terms  of  the 
local  position  but  also  the  time-step  size  At,  such  that 


xw  = XwU,?,At)  = x(At)/A^U,?) 

X0  = X0(5,S,At)  = x( At)/ A^( e , c)  • 


(6.9) 


CHAPTER  VII 
CONCLUDING  REMARKS 


The  adaptive  grid  generation  based  on  the  variational 
principle  was  investigated  in  detail  for  the  projectile 
aerodynamics  problem.  The  one-dimensional  principle  has 
been  successfully  applied  to  adaptive  boundary  grid 
generation,  as  well  as  the  two-dimensional  variational 
principle  to  the  grid  generation  in  the  domain.  Excellent 
consistency  has  been  achieved  with  the  same  corresponding 
penalty  parameters.  The  developed  one-dimensional  grid 
generator  can  be  applied  to  any  complex  configuration  of 
boundary  such  as  a curvilinear  surface  combined  with  a 
polygonal  shape. 

The  penalty  parameters  for  orthogonality  and 
concentration  of  grid,  treated  as  constants  in  the  work 
reported  to  date,  have  been  investigated  to  show  that  they 
are  of  the- same  order  of  magnitude  and  can  be  treated  as 
variable.  These  variables  are  in  terms  of  not  only  the  cell 
grid  spacing  over  the  entire  domain  but  also  the  time-step 
size.  The  gradient  of  pressure  is  chosen  as  the  control 
weight  function  in  the  computation  of  inviscid  transonic 
flow  problem.  However,  the  generated  grid  spacing  in  the 
normal  streamwise  direction  close  to  the  body  surface  is 
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often  too  large  and  causes  the  numerical  instability  of 
solution  algorithm.  In  order  to  solve  this  problem,  the 
investigation  has  shown  that  an  artificial  clustering  of  the 
grid  points  close  to  the  surface  boundary  can  be  used  to 
ensure  the  convergent  and  stable  computation  in  the  Navier- 
Stokes  code. 

The  adaptive  grid  generation  code  coupled  with  an 
axisymmetric  version  of  the  thin-layer  Navier-Stokes  code 
has  been  used  successfully  for  the  computation  of  the 
transonic  projectile  flow  field.  The  coupled  computer  code 
offers  an  automatic  calculation  procedure  for  changing  the 
time  step,  penalty  parameters,  and  the  adaptive  grid 
network. 

Comparisons  were  made  of  the  results  of  computations 
from  adaptive  and  fixed  grid  systems  with  experimental  data 
for  projectile  inviscid  transonic  flow  problem.  The 
comparisons  show  that  the  adaptive  grid  system  gives  a 
better  agreement  with  the  experimental  data  than  the  fixed 
grid  system. 


APPENDIX  A 

NOMENCLATURE  OF  INPUT  PARAMETERS 
IN  THIN-LAYER  CODE 

A . 1 NMAX, 

JMAX,  KMAX,  LMAX,  ND 

NMAX: 

The  maximum  number  of  time  steps 

JMAX: 

The  maximum  number  of  grid  points  along  5 
direction 

KMAX: 

The  maximum  number  of  grid  points  along  j) 
direction 

LMAX : 

The  maximum  number  of  grid  points  along  5 
direction 

ND: 

The  number  usually  equal  to  KMAX 

A. 2 DT,  DX,  DY,  DZ,  FSMACH,  SMU,  ALP 


DT: 

Time  step  At 

DX: 

Step  size  Ax  of  x-direction 

DY: 

Step  size  Ay  of  y-direction 

DZ: 

Step  size  Az  of  z-direction 

FSMACH:  Free  stream  Mach  number  1M 


SMU: 

Explicit  smooth  parameter 

ALP: 

Angle  of  attack  a 
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A. 3 IREAD, 

I WRIT , NGRI,  INVISC,  NP 

IREAD: 

Number  of  readings 

0:  starting  from  N=1  to  N=NMAX 

1 : restarting  from  N=NMAX+1 

IWRITE: 

Number  of  writings  to  control  output 

NGRI : 

Number  of  grid  input  data  types 
1 : 2-dimensional 

2:  Tape  3 

INVISC: 

Number  of  flow  types 
0:  Inviscid  flow 

1 : Viscous  flow 

NP: 

Number  of  time  step  intervals  to  print  out  the 
output 

A. 4 CNBR,  OMEGA,  PDOV 


CNBR: 

Courrant  number 

OMEGA: 

Angular  speed  w 

PDOV: 

A factor  transferring  the  Mach  number  to  spin 
number  (i.e.,  SPIN=PDOV*FSMACH) 

A. 5 RE,  PR,  RMUE,  TZ 


RE: 

Reynold's  number  Re 

PR: 

Prandtl  number  Pr 

RMUE: 

Coefficient  of  viscosity 

TZ: 

Stagnation  temperature  T0 
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A. 6 LAMIN,  RM 

LAMIN:  Number  of  controls  of  flow  case 

0:  Laminar  flow 

1 : Turbulent  flow 

RM:  Implicit  smooth  parameter  ej 


APPENDIX  B 

PARTIAL  PROGRAM  LISTING 


The  grid  generation  code  "GRIDGEN"  can  be  implemented 
by  adding  several  subprograms  so  that  the  user  will  have  an 
option  to  generate  an  adaptive  grid  network  for  the  accurate 
solutions.  The  implementation  and  modification  are 
described  in  this  Appendix.  In  the  process,  efforts  have 
been  made  to  keep  the  changes  as  minimal  as  possible.  In 
fact,  only  seven  subroutines  ( ADGRD , GRADP , WFCN1 , ADBD , WFCN2 , 
DIFF,  and  RELAX)  with  executable  statements  are  added  to 
GRIDGEN  to  generate  adaptive  grid  points  on  the  boundary  and 
in  the  domain.  The  modifications  made  to  GRIDGEN  and  the 
relevant  information  required  for  using  adaptive  GRIDGEN  are 
described  in  the  following: 

A.  Input  Data 

Two  data  cards  are  added  at  the  beginning  of  the  input 
data  file  of  GRIDGEN,  and  their  executable  statements  now 
read  these  two  additional  input  data  cards  as  Hsu  [8]  used 
the  statements  in  the  program  MAIN  of  GRIDGEN.  The  first 
datum  read  is  IADAPT,  such  that 


r 0 


IADAPT  < 


V 1 


implemented  by  GRIDGEN 
implemented  by  adaptive  GRIDGEN 
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The  second  data  read  are 

IE(1),IE(2),IE(3)»IE(4)»IE(5)»IE(6),IE(7)»IE(8) 

These  data  provide  information  and  instruction  for  the 
change  of  number  of  grid  points  on  each  segment  of  inner 
boundary  (i.e.,  the  change  of  boundary  grid  spacing).  For 
example,  there  are  23  points  on  ogive,  22  points  on  the 
cylinder,  17  points  on  the  boattail,  and  8 points  on  the 
sting.  We  then  set  IE ( 1 ) =1 , IE(2)=23,  IE(3)=23,  IE(4)=45, 
IE(5)=45,  IE(6)=62,  IE(7)=62,  and  IE(8)=70.  The  use  of  the 
COMMON  IADAPT  and  COMMON  IE(8)  are  available  for  the 
relevant  subprograms  such  as  they  are  used  in  ADGRD  to 
generate  the  adaptive  grid  points  on  each  segment  of 
boundary.  The  program  then  continues  as  before  to  generate 
inner  boundary  grid  points  by  calling  appropriate 
subroutines. 

B.  Inserted  Statements 

In  order  to  add  a new  option  for  adaptive  grid 
generation  into  GRIDGEN.  It  is  suitable  to  insert  two 
statements  in  the  subroutine  "ALGRD,"  the  first  one  is  the 
"IF"  statement  for  the  optional  number  IADAPT=1 , then  we  go 
to  the  second  statement  "CALL  ADGRD,"  otherwise  we  bypass 
this  executable  statement.  Since  we  can  use  directly  the 
implementation  of  clustering  grid  points  close  to  the 
boundary  after  the  adaptive  grid  generated  in  the  subroutine 
ALGRD  of  GRIDGEN,  this  is  why  we  insert  the  execution 
statement  for  adaptive  grid  generation  initially  in  ALGRD. 
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C.  Inserted  Subprograms 
ADGRD : 

This  subroutine  is  an  adaptive  grid  solver  called 
when  IADAPT=1 , the  subroutines  GRADP,  ADBD,  and  RELAX 
are  called. 

GRADP: 

The  calculations  of  pressure  gradients  or  any  other 
gradients  with  one  dimension  only  for  the  boundary  will 
be  executed  in  preparation  for  constructing  the  weight 
function.  Also,  the  calculated  pressure  gradients  are 
modified  by  applying  Eqn.  (5-4). 

WFCN1  ( ) : 

The  one-dimensional  boundary  weight  function  is 
executed  by  applying  Eqn.  (5.2)  with  the  modifying  Eqn. 
(5.6) . 

ADBD( ) : 

This  subroutine  is  the  calculation  of  adaptive 
boundary  grid  points  executed  by  Eqn.  (2.64).  Prior  to 
implementation  of  the  Eqn.  (2.64),  the  Cartesian 
coordinates  are  transformed  into  curvilinear 
coordinates,  thus,  the  boundary  points  may  be  described 
in  one  dimension.  We  then  apply  the  Newton-Raphson 
method  with  SOR  to  calculate  the  adaptive  grid  on  the 
boundary.  Also,  we  have  to  transfer  the  curvilinear 
system  back  to  the  Cartesian  system  by  employing  cubic 
spline  function  as  the  interpolating  function  to  find 
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new  (x,z)  values  on  the  boundary.  This  subroutine  is 
called  by  ADGRD  five  times  (i.e.,  an  outer  boundary  plus 
four  segments  combined  to  an  inner  boundary  of 
projectile  surface).  The  subroutine  WFCN1  is  called. 
WFCN2(-) : 

The  two-dimensional  weight  function  is  calculated 
by  Eqn.  (5.1)  with  the  modifying  Eqn.  (5-4)  and  (5.5). 
Also,  the  Jacobian  determinant  is  calculated  after  each 
iteration  in  RELAX.  The  derivatives  of  weight  function, 
Wx  and  Wz,  with  transformations  are  also  calculated. 
DIFF( ) : 

This  subprogram  calculates  the  metrics,  second 
derivatives  Jacobian  determinant,  coefficients  of 
highest  derivative  terms  and  residuals  of  the  coupled 
Eqn.  (2.36) . 

RELAX  ( ) : 

This  subprogram  is  the  solver  for  two-dimensional 
grid  calculations  executed  by  the  N-R  method  with  SOR  in 
using  Eqn.  (2.57)  and  (2.58).  The  subroutine  WFCN2  is 
called  for  each  iteration  since  the  Jacobian  determinant 
should  be  changed  currently.  The  subroutine  DIFF  is 
called  to  calculate  either  the  adaptive  grid  points  if 
IADAPT=1  or  regular  grid  points  if  IADAPT=0.  The  last 
task  is  to  extrapolate  linearly  the  interior  grid  lines 
orthogonal  to  the  axis  and  downstream  boundary.  The 
iteration  continues  until  it  converges  up  to  either  the 
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allowable  control  residual  RES=0.01  or  the  specified 
maximum  iteration  number  ITERM=500.  Note  that  this  new 
subroutine  RELAX  is  replaced  by  the  old  "RELAX"  in 
GRIDGEN. 


D.  The  Flow  Program 

The  summarization  of  above  descriptions  results  in  the 
flow  program  for  adaptive  grid  calculation  as  shown  in  Fig. 

A.  1 . 


Return 


Fig.  A.1  Flow  program  for  adaptive  grid  calculation 


E.  Computer  Subprograms 

The  computer  subprograms  are  shown  in  the  following 


pages. 


O O O O O O O O O o o 
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SUBROUTINE  AOGRQ 

IMPLICIT  00U3L5  PRECISION  *6(A-H,0-Z) 

DOUBLE  PRECISION  *6  L30/L3V 

COMMON  JMAX,  K MA  X , JM,  KM/  N30D/  J30D 

COMMON/ GRD/X1 ( 90, 100 > , Y1 (90, 100) 

COMMON/ VARS/ 0(90,6/100) 

COMMON/VAR3/P(60/100) 

^ COMMON^ / SO  UOY/^XXU  00),  YY(IOO)/  XSdOO),  YS(IOO),  SS(IOO),  S(100) 
COMMON  /OUTER/  JAC8),  JS(8) 

rnMuo!!  ',ZI  CA0<8>'  CA1C8),  CA2  (8 ) / CBO  (8)  / C31C8),  C32(3),  CARC(S) 
COMMON  /ARRAY/  AC100),  3(100),  C(100),  0(100),  F(100),  HdOO) 
COMMON  /FIX/  XFIX,  YFIX 

COMMON  /ELLIP/  R X XO , R X YO , R Y X D , R Y YO , R X , R Y , R D ET , LB  0, L 3 V , L 3 S 
COMMON/PYHSU/I ADAPT 
COMMON/ C0UNT2/NN/ISCAL 
COMMON  /PARA/  AFA/STA 

COMMON/ WEIGHT/ W ( 90, 1 00), DWX (90,1 00), DWY (90,  100) 

CCMMON/SEGP/IE (8) 

COMMON /GGONLY/IGG/NGG/NGM AX ,MM 
AFA=1 . E-6 
3TA=1 .0 
EPS=1 . D-4 


CALL  OUTPUT  TO  FIND  THE  PRECEDING  APPROXIMATE  SOLUTION 
CALL  0UTPUT(0/0/0/ISCAL/1) 

FIND  THE  PRESSURE  GRADIENTS  ALONG  THE  30UN0ARY 
CALL  GRADP 

CALCULATE  THE  ADAPTIVE  GRID  ON  INNER  BOUNDARY 

NOTE  THAT  THERE  ARE  FOUR  SEGMENTS  CCMSINED  THUS  WE  MUST 

CALCULATE  SEPARATELY 

N1=IE(1 )+1 

N2=IE(2)-1 

CALL  ADBDdTERM,  EPS, WM EGA, N1,N2, 1,1, 1,1,1) 

DO  43  J=N1,N2 
XI  ( J,1 )=T(J) 

Y1 ( J, 1 ) =TS ( J ) 

43  CONTINUE 
N3=IE (3 ) +1 
N4=IE(4)-1 

CALL  ADBD( IT ER M, EPS, WMEGA/N3/N4, 1,1, 1,1,1) 

DO  44  J=N3,N4 

44  XI (J,1 ) *T ( J ) 

N5=IS(5)+1 

N6«IS(6)-1 

CALL  ADBDdTERM, EPS, WMEGA/N 3, N6, 1,1, 1,1,1) 

DO  46  J=N5,N6 
XI  (J,1)=T(J) 

Y1  (J,1 )=TS(J) 

46  CONTINUE 


o o o o o o 
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N7=IE  C7)+1 

CALL  A03DCITERM, EPS, WMEGA,N7,JM, 1,1, 1,1,1) 

DO  48  J=N7,JM 
XI (J,1 )=T(J) 

48  CONTINUE 

CALCULATE  THE  ADAPTIVE  GRID  ON  OUTER  BOUNDARY 

N J = J B ( 1 ) 

NP=N J+1 
NM=NJ-1 
5PM=1 .D-2 

CALL  ADB0(ITERM,EPM,WMEGA,2,JM,1,KMAX,KMAX,1,1) 
DO  54  J=2,JM 
XI (J,KMAX)=T(J) 

Y 1 (J,KMAX)=TS(J) 

54  CONTINUE 

CALCULATE  THE  ADAPTIVE  GRID  IN  THE  DOMAIN 


RE3D=1 . D-2 

CALL  RELAX(ITERM,RESD,WMEGA,2,JM,1,2,KM,1) 

WRITE  (16,200)  C(X1(J,K),J=1,JMAX),K=1,KMAX),((Y1 (J,K),J=1,JMAX), 
1 K=1,KMAX) 

IFCMM  .GE.  2)  NCLUS=3 
RETURN 

200  F0RMAT.(14F9.5) 

END 


o o o o o o 
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SUBROUTINE  OR A OP 

IMPLICIT  DOUBLE  PRECISION  *6(A-H,0-Z) 

COMMON  JMAX/  K MA  X / JM/  KM/  N300/  J30D 
C0MM0N/VAR3/P (60/100) 

COMMON/ V A RS/Q (90/6/1 00) 

COMMON  /GRAD/  OPDX(60/100)/OPDY(60/100) 
COMMON/GANDM/RGAMA/RMACH 

CALCULATE  THE  PRESSURE  GRADIENTS 

DO  20  K=1/K MAX/KM 
00  10  J=2/JM 

0P0X(K/J)=0.5*(P(K/J+1)-P(K/J-1)) 

10  CONTINUE 

DPDX(K/1)=0.5*(-3.*P(K/1)+4.*P(X/2)-P(K/3)) 

DPDX(K/JMAX)=0.5*(3.*P(K/JMAX)-4.*P(K/JM)+P(K/JM“1)) 

20  CONTINUE 

TAKE  AVERAGE  AT  THE  CHANGED  SIGN  OP  PRESSURE  GRADIENTS 

DO  30  K = 1 /KMAX/KM 
DO  30  J=2/JM 
JP=J+1 
JR=J-1 

PXR=P(K/J)-P(K/JR) 

PXP=P(K/JP)-P(K/J) 

XPR=PXP*PXR 

'0  CONTINUE1*1"*  DPOX  (K/J)=0.5*(DPCX  (K/JP)-DPDX(K/JR) ) 

RETURN 

END 


° ° n n o o o 
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SUBROUTINE  WFCN1(ITER/J1/J2/JV,K1/K2/KV/JO) 

IMPLICIT  DOUBLE  PRECISION  *6(A-H,0-Z) 

COMMON  J M A X / KMAX/  JM/  KM/  N300/  JBOD 
COMMON/GRO /XI (90/ 100), Y1 (90,100) 

C0MM0N/VAR3/?(6O/10O) 

COMMON / VA RS/ C ( 90/6/100) 

COMMON/ W EIGHT/ W ( 90/ 1 00) / DWX ( 90/1 00 >/OWY( 90/1 00) 

COMMON  /ELLIP/  R X X 0 / R X Y D / R Y XD, R Y Y 0 , R X / R Y / R 0 E T, L 3 0/ L B V , 15 S 
COMMON  /ARRAY/  A(100),  5 (1  00)/  C(1  00)/  0(100)/  F(100)/  HdOO) 
COMMON  /BOUDY/  XX(IOO)/  YYdOO),  XS(IOO)/  YS(IOO)/  SS(IOO)/  SdOO) 
1 / T (1 00) / TS (1 00) 

COMMON  /GRAD/  DPDX (60,1 00), OPDY (60,1 00) 

COMMON  /PARA/  AFA/5TA 

I M A X = 3 

K = K 1 

JI=J1 -1 

JF=J2+1 

TRANSFER  TO  CURVILINEAR  SYSTEM  FOR  ONE-DIMENSION  ON  BOUNDARY 

IF ( ITER  .GT.  1 ) GO  TO  40 
SS( JI)=0.0 
DO  40  J=J1,JF,JV 
JR= J-1 

SS(J)  = SS(JR) +DSCRT  ( (XI  ( J/K) -XI  (JR,K))**2+(Y1  (J/K)-Y1  ( JR,K) ) **2) 

40  CONTINUE 

TRANSFER  THE  PRESSURE  GRADIENT  TO  WEIGHT  FUNCTION 

DO  60  J=J1/J2/JV 

SXD  = 0. 5 * ( S 3 (J  + 1 )-SS  (J-1 ) ) 

W(K/J)=OPDX(K/J) /SXD 
W(K,J)=AFA+5TA*DA3S(W(K,J)) 

60  CONTINUE 

SXD1=0.5*(-3.*SS(JI)+4.*S3(J1)-SS(J1+1)) 

SXDM=0.5*(3.*SS(Jc)-4.*3S(J2)+SS(J2-1)) 

W(K,  JI)=DPDX  (K/JD/SXD1 
W(X/ JF) =DPDX (K/JF ) /SXOM 
W(K,JI)=AFA+5TA*DA5S(W(K,JI)) 

W(X/JF)=AFA+3TA*0ABS(W(K/JF)) 

CALCULATE  THE  WEIGHT  cUNCTION  SMOOTHLY 

IF (X  . EQ.  KMAX ) I M A X = 2 0 

DO  55  I=1 /I MAX 

DO  80  J=J1/J2/JV 

WW  = 0. 5* (W(K/ J + 1 ) +W(K/ J-1  ) ) 

W(K,J)=W(iK/J)+0.4*(WW-W(K/J)) 

80  CONTINUE 
85  CONTINUE 

CALCULATE  THE  DERIVATIVE  OF  WEIGHT  FUNCTION 
DO  90  J=J1 /J2/JV 

DWX(K/J)=0.5*(W(K/J+1)-w(K/J-1)) 

90  CONTINUE 
RETURN 
END 
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SUBROUTINE  AD3D(ITERM,RES, WMEGA, J1, J2, JV,K1 ,K2,KV, JO) 

IMPLICIT  DOUBLE  PRECISION  *6(A-H,0-Z) 

DOUBLE  PRECISION  *6  J C 03 / L A X, L 3 Y , L 3 X , LC Y , LB  0/ L3 V / LB S 
CO’MMON  JMAX/  KMAX,  JM,  KM/  NSOO/  J50D 
COMMON/ GRO /XI (90,100),Y1 (90/100) 

COMM ON /VARS/QC 90/6/100) 

C0MM0N/VAR3/PC60/100) 

COMMON/ WE IGHT/ W( 90/1 00)/ DWX (90, 100), OWY (90/100) 

COMMON  /ELLIP/  R X XD , R X YO , R Y X 0/ R Y Y 0 / R X, R Y, RD ET / LB  0 , LS V , L 5 S 
COMMON  /ARRAY/  A(100)y  3(100),  C(ICO),  0(100),  F(100),  H(100) 
COMMON  / 30U0 Y/  XX(IOO),  YY(IOO),  XS(IOO),  YSdOO),  SSdOO),  SdOO) 
1 / T ( 1 00 ) / T S (100) 

COMMON  / D A R A / A F A , 3 T A 

CALCULATION  FOR  ADAPTIVE  BOUNDARY  GRID  POINTS  3Y  USING 
NEWTON-R APHSON  ITERATION  WITH  SUCCESSIVE  OVER  RELAXATION  METHOD 


BEGIN  TO  ITERATE 
DO  60  ITER=1 /ITERM 
ITCON  = 1 

CALL  WFCN1 (ITER/J1/J2/JV/K1 ,K2/KV,J0) 

K = K 1 
JI= J 1 - 1 
JF=J2+1 

NORMALIZE  THE  LUM3DA  W BY  DIMENSION  ANALYSIS 
WT=Q.O 

DO  20  J = J1*JF 

20  WT=WT+0.5*(W(K/J)+W(K/J-1))*(S3(J)-SS(J-1)) 

TP=JF-JI 
TL=SS (JF) 

SL=(TL/TP)**3 
WBAR=WT/TL 
SL3V=W3 AR*SL 
R3V=L3V/ SL3V 

I F ( R B V .GT.  1.06)  RBV  = 1.D6 

CALCULATE  THE  ADAPTIVE  30UNDARY  POINTS 

DO  40  J = J 1 / J2, JV 
JP= J+1 
J R=  J - 1 

SXD  = 0.5*(SS( JP)-SS(JR)  ) 

SXXD=SS(JP)-2.*SS(J)+SS(JR) 

RS=SXXD+R3V*(SXD**4)*DWX (K,J) / (2. * ( 1 . + R S V * ( S XD** 3 ) *W ( K, J ) ) ) 
RSXD=-2.0 
RRS=RS/RSXD 
SS(J)=SS(J)-WMEGA*RRS 
C CONTROL  THE  ACCURACY  OF  SOLUTION 

ARS=DA3S(RS) 

I F ( A R S .GT.  RES)  ITCON  = 0 
40  CONTINUE 

IF(ITCON.NS.O)  GO  TO  S3 
60  CONTINUE 
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80  CONTINUE 

TRANSFER  3 ACK  TO  CARTESIAN  SYSTEM  3 Y USING  CUBIC  SPLINC 
INTERPOLATION 

IFCK  . = 2.  KMAX ) GO  TO  120 

CALCULATE  THE  X/Y  VALUES  ON  INNER  BOUNDARY 

S(JI)=0.0 
DO  100  J=J1/JF 

100  S(J)=S(J~1)+DS0RT ( (XX(J)-XX (J-1) )**2+(YY(J)-YY(J-1 ) )**2) 
^ ^ ^ ^ CSPLIN(SS/TrS/XX/A/B/C/D/FrHrJ1/J2/’JI/JF) 

*^^CL  0SPLIN(SS/TS/S/YY/A/3rC/D/F/H/J1/J2^JI/JF) 

GO  TO  160 

CALCULATE  THE  X/-Y  VALUES  ON  OUTER  BOUNDARY 

120  CONTINUE 
S(JI)=0.0 
DO  140  J = J 1 / J F 
JR=J-1 

S(J)=S(JR)+DSQRT((XS(J)-XS(JR))**2+(YS(J)-YS(JR))**2) 

140  IF(SCJ)  . E3.  S ( J R ) ) S ( JR ) =0 . 5 * ( S ( J ) + 3 C J R- 1 ) ) 

CALL  CSPLIN(SS/T/S/XS/A/3/C,D/P/H/J1,J2/JI/JF) 

^ ^ 4 L CSPLINCSS,TS/S,YS,A,3,C,D/F/H,J1,J2,JI,JF) 

160  CONTINUE 
TL=SS(JF) 

WRITE  (6r  200)  ITER  r TL 

WRITE (6/ 1 90)  <WCK,J>/J«JI,JF> 

190  FORMAT (12E11.4) 

200  F0RMATC//5X/ "THE  NUM3ER  OF  ITERATIONS  = "/T4,6X/ 

1 "TOTAL  LENGTH  = ",=1 1 .4/) 

RETURN 

END 


o o o n 


137 


SUBROUTINE  RELAX ( ITER M/ RES/ WMEGA,J1,J2/JV,K1/K2/KV) 

IMPLICIT  DOUBLE  PRECISION  *6(A-H/0-Z) 

DOUBLE  PRECISION  *6  JC03/LAX/L3Y/L3X/LCY/L30/L3V/L3S 
COMMON  J M A X / K M A X / JM/  KM,  N300,  J30D 
COMMON/GRD/Xl (90/100)/Y1 (9C/100) 

COMMON/ VARS/ 3 (90/ 6/1 00) 

COMMON/ V4R3/=( 60/ 100) 

CO MM ON/ W EIGHT/ W(90/10Q)/DWX(90/10Q)/DWY (90/100) 

COMMON  /3L0CK0/  XXF(S0/40)/XEF(30/40)/YXF(30/40),YEF(80/40) 

COMMON  /ELLI?/  RXXD/RXY0/RYXD/RYYD/RX/RY/RDET,L30/LBV/LaS 

COMMON/PYHSU/IAOAPT 

C0MM0N/C0UNT2/NN/ISCAL 

COMMON  /ARRAY/  A(IOO)/  3(100)/  C(1C0)/  0(100)/  F(100)/  H ( 1 n Q ) 
COMMON  / 30UD  Y/  XX(IOQ),  YY  (100)/  XSdOO),  YSdCO)/  SSdOO)/  SdOO) 
1 / T (1 00) / TS (100) 

COMMON  /PARA/  AFA/3TA 
C 

NEWTON-RAPHSON  ITERATION  WITH  SUCCESSIVE  OVER  RELAXATION  METHOD 

X0=X1(1,1) 

JI=J1-1 

JF=J2+1 

KI=K1-1 

KF=K2+1 

IF (IAOAPT  . EQ.  0)  R30=1.0' 

IF(IAOA?T  .33.  0)  R 3 V = 0 . 0 
C BEGIN  TO  ITERATE  BY  USING  ADI  METHOD 

DO  90  I T E R = 1 /ITERM 
I T C 0 N = 1 

IF (I  ADAPT  .33.  0)  GO  TO  10 
CALL  WFCN2  (ITER) 

10  CONTINUE 

IF(MC3(ITER/2)  .33.  0)  GO  TO  45 
DO  40  K=K1,K2/KV 

ITERATING  CALCULATION  ALONG  THE  XI-LINE 
DO  40  J=J1/J2,JV 

NORMALIZE  THE  LUMEDA  V AND  LUM3DA  0 BY  DIM=NSI0N  ANALYSIS 
IF  (IAOAPT  .33.  0)  GO  TO  35 
0SX=DS3RT((XXF(J,K))**2+(YXF(J,K))**2) 

S L = D S X * *4 
SL3V=W (K/J) *SL 
R3V=L3V/SLBV 
R30=L30/SL 
35  CONTINUE 

CALL  DI=F(J,K/RE0/R3V) 

RRX=(RX*RYYD-RY*RXYC)/RDET 
RRY=(RY*RXXD-RX»RYXD)/ROET 
XI (J,K) =X1 (J/K)+WMEGA*RRX 
Y1(J/K)=Y1(J,K)+WMEGA*RRY 
CONTROL  THE  ACCURACY  0=  SOLUTION 
ARX=DA3S(RX) 

ARY  = D A3S  (R Y) 
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IFCARX.GT .RES. AND. ARY. GT .RES)  ITC0N=0 
40  CONTINUE 
GO  TO  70 
45  CONTINUE 

DO  70  J=  J1  / J 2/  JV 

ITERATING  CALCULATION  ALONG  THE  ETA-LINE 
DO  70  X=X1,X2,XV 

NORMALIZE  THE  LUM3DA  V ANO  LUM5DA  0 3Y  DIMENSION 
IF (IAOAPT  . EQ . C)  GO  TO  60 
D3E=DSgR7((XSF(J,X))**2+(YEF(J,X))**2) 

SL=DSE**4 
SL3  V = W ( X , J ) * S L 
R3V=L3V/SL3V 
R30=L30/3L 
60  CONTINUE 

CALL  DIFF(J,X,R5C,R3V) 

RRX=(RX*RYY0-RY*RXYD)/ROET 
RRY=(RY*RXXO-RX*RYXD)/R0ET 
X1(J,X)=X1  (J,X)+WMEGA*RRX 
Y 1 (J,X)=Y1  (J,X)+WMSGA*RRY 
CONTROL  THE  ACCURACY  0=  SOLUTION 
ARX=A3SCRX) 

A R Y = A 3 3 CRY) 

IFCARX.GT.RE3.AND.ARY.GT.RES)  ITC0N=0 
70  CONTINUE 

IFCITCON.NE. 0)  GO  TO  100 
X?  = 1 

XT  = X1 (2,2) 

IFCXT  .LT.  XO)  GO  TO  35 
75  XPsXe-M 

XT=X1 ( 2, XP ) 

IFCXT  .LT.  XO)  GO  TO  SC 
GO  TO  75 
30  CONTINUE 

7X1 = X1 C2,XP) -XO 
T X 2 = X 1 ( 2 , X P ) -X 1 C 2 , 1 ) 

DC  35  K = 2, K P 

0X2  = X1  (2/O-XI  (2 , X-1  ) 

XI  (1,X)=X1  (1,K-1)+TX1*(DX2/TX2) 

Y1 C JMAX,K)=Y1 ( JM, X) 

25  CONTINUE 
XP=X?*1 
DO  90  X = X P , .< M 
XI  ( 1 , X ) = XI ( 2 , X ) 

Y 1 ( JMAX,K)=Y1  C J M , X ) 

90  CONTINUE 
100  CONTINUE 

•JRIT  = (6,  200)  ITcR,L30,L3V 

200  FORMAT (/ /5X,"THE  NUM3ER  0 = ITERATIONS  = ",I4/6X, 

1 "LU.M5  0 = " , E 1 1 , 4, 5 X , "LUM3  V = ",=11.  4//)* 

RETURN 

END 
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SUBROUTINE  31 FF ( J ,K / R SO, R 3 V ) 

IMPLICIT  DOUBLE  PRECISION  *6(A-H,0“Z) 

DCU3LE  PRECISION  *6  JC03 / LA  X/ L3 Y / L 3X / LC Y/ L30/ L3 V / LBS 
COMMON  UMAX/  UMAX/  JM/  KM/  N3  0D/  J300 
COMMON/ GRD/X1 (90/100) /Y1 (90/100) 

COMMON/ VAR S/Q( 90/ 6/ 100) 

C0MM0N/VAR3/PC60/100) 

COMMON/WEIGHT/W(9U/100)/DWX(90/100)/DWY(9 0/100) 

COMMON  / BLOCKO / XXF(80/40)/XEF(30/40)/YXF(30/40)/YEF(80/40) 

COMMON  /ELLIP/  R X X D / RX Y 0 / R Y XD / R Y Y 0 / R X / R Y / RD ET/ L 3 0/ l* V/ L3 S 

COMMON/PYHSU/IADAPT 

COMMON  /PARA/  AFA/BTA 

JP=J+1 

JR=J-1 

KP=K+1 

KR=K-1 


CALCULATE  THE  !=IRST  DERIVATIVE  TERMS 

XXD=(X1 (JP/K)-X1 (JR,K))/2.C 
X ED= ( XI (J/KP)-X1 (J/KR))/2.0 
YXD=(Y1 (JP/K)-Y1 (JR/K) )/2.0 
YED= ( Y1 (J/KP)-Y1 (J,KR))/2.C 

CALCULATE  THE  SECOND  DERIVATIVE  TERMS 

XXXD=X1 (JP/K)-2.*X1 (J/K)+X1 ( JR/K) 
XESD=X1CJ/KP)-2.*X1CJ/K)+X1(J/KR) 

YXXD  = Y1  (JP/K)-2.*Y1(J/K)+Y1(JR/K) 

YEED=Y1 (J/KP)“2.*Y1(J/K)+Y1 (J/KR) 

XXED=(Xl(JP,KP)-X1(JP/KR)-X1(JR/KP)+X1(JR/KR))/4.0 

YXED=(Y1(JP/KP)-Y1(JP/KR)-Y1(JR/KP)+Y1(JR/KR))/4.0 

CALCULATE  THE  JACOBIAN  DETERMINANT 

JCOE-XXD*YEO”XED*YXD 
I F C JC0B.EQ.0.0)  GO  TO  140 

CALCULATE  THE  VALUES  OF  SOURCE  TERMS  OF  THE  EQUATION 
THE  DECISION  OF  THE  ADAPTIVE  GRID  GENERATION 

IFCIADAPT  .EQ.  0)  GO  TO  40 
S0X=0,5*JC03*DWX(K/J) 

S0Y=Q.3*JC03*DWYCK/J) 

40  CONTINUE 

CALCULATE  THE  COEFFICIENTS  OF  HIGHEST  DERIVATIVE  TERMS 

ALFA  = XED**2  + YED**  2 
BETA=-2.*(XXD*XED+YXD*YED) 

GAMA=XXD**2+YXD**2 

AA=(YXD**2+YED**2)/(JC03**3) 

33=-(XXD*YXD+XED*YED)/(JC03**3) 

CC=(XXD**2+XED**2)/ (JC03**3) 
CAA=AA*ALrA+R30*XED*XED+R3V*YED*YED*W(K/J) 
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CA3  = AA*3ETA  + R30*2.*C2.*XX0*XED  + YX3*YED)-R3V*(2.*YX0*YED)*W(K,rJ) 

CAC=AA*GAMA+R30*XX3*XXD+RSV*YXD*YXD*W(K,J) 

C3A=33*ALFA+R30*XED*YED-RBV*XED*YED*W(K/J) 

CB3=33*BSTA+R30*CXX0*YED+XED*YXD)+RSV*CXXD*YED+XE0*YXD)*W(K,J) 

C3C=33*GAMA+R30*XX3*YX0-R3V*XXD*YXD*W(K,J) 

CCA=CC*ALFA+R30*YED*YED+R3V*XE0*XEC*W(K/J) 

CC3=CC*3STA+R30*2.*(XXD*XED+2.*YXD*YED)-R3V*C2.*XXD*XED)*WCK,J) 

CCC=CC*GAMA+R30*YXD*YXD+R3V*XXD*XXC*WCK,J) 

LAX=CAA*XXXD+CA3*XXED+CAC*XEED 

L3Y=C3A*YXXD+C33*YXED+C3C*YEED 

L3X=C3A*XXXD+C33*XXED+C3C*XEED 

LCY=CCA*YXXD+CCB*YXED+CCC*YEED 

CALCULATE  THE  RESIDUALS  OF  THE  COUPLED  EQUATIONS 

RX=LAX+L3Y+RSV*S0X 

RY=L3X+LCY+R3V*S0Y 

RXXD=-2.*(CAA+CAC) 

RXYD=-2.*CC3A+C3C) 

RYXD=~2.*CC3A+C3C) 

RYYD=-2.*(CCA+CCC) 

RDET=RXYD*RYXD-RXXD*RYYD 
120  CONTINUE 
RETURN 

140  WRITE  ( 6/  1 30) 

STOP 

130  FORMAT (/ /5X/" JAC03I  DETERMINANT  IS  SINGULAR") 

END 
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SU3R0UTINE  WFCN2  (ITER) 

IMPLICIT  00U3LE  PRECISION  *6(A-H/C-Z) 

COMMON  UMAX/  K MA  X / JM/  KM/  N300/  J30D 

C0MM0N/GRD/X1(90/100)/Y1 (90/100) 

C0MM0N/VAR3/P(60/100) 

COMMON/ VARS/ 3 (90/ 6/ 100) 

COMMON/WEIGHT/ W (90/1 00 )/DWX(90/100)/DWY(90/100) 

COMMON  / 3L0CKD/  X X F ( 30/4 0 ) / XEF ( 30 / 40 ) / YXF ( 30/ 40 ) / YEF ( 3 0/40 ) 

COMMON  /ELLIP/  R XXO/ RXYD / R Y XO/ R Y YD/ RX/ R Y/ RDET/ L3 0/ L3 V/ L3S 
COMMON  /ARRAY/  A(100)/  3 (100)/  C(1C0)/  0(100)/  FdOO)/  HdOO) 
COMMON  / 30U0  Y/  XX(IOO)/  Y Y (1  00)  / XSdOC)/  YS(IOO)/  SS(IOO)/  S(100) 
1 / T d 00 ) / TSdOO) 

COMMON  /GRAD/  DP 0 X ( 6 0/ 1 0 0 ) / 0 PD Y ( 6 0 / 1 00 ) 

COMMON  /PARA/  AFA/3TA 
COMMON /GANDM/R GAMA/ RMACH 
IMAX=3 

CALCULATE  THE  METRICS  AND  JAC03IAN  AT  EACH  ITERATION  FROM  RELAX 


KV=1 

DO  20  K = 1 /KMAX/KV 
DO  10  J=2/JM 

XXF(J/K)=(X1 (J+1/K)-X1 ( J-1 / K) ) /2 . 0 
YXF(J,K)=(Y1(J+1,K)-Y1(J-1/K))/2.0 
10  CONTINUE 

XXF(1/K)=(-3.*X1 (1/K)+4.*X1 (2/K)-X1(3/K))/2.0 
YXF(1/K)=(-3.*Y1(1/K)+4.*Y1(2/K)-Y1(3/K))/2.0 
XXF(JMAX/K)  = (3.*Xl(JMAX/K)-4.*X1(JM/K)t>X1(JM-1/K))/2.0 
YXr ( JMAX/K) =(3 .*Y1 ( JMAX/K)-4.*Y1 ( JM/ K)+Y1 ( JM-1 /K) ) / 2.0 

20  CONTINUE 
JV  = 1 

DO  40  J = 1 /JMAX/JV 
DO  30  K = 2 / K M 

XEF(J/K)  = (X1  (J/K  + D-  X1  (J/K-1  ) ) / 2.0 
YEF(J/K)=(Y1 (J/K+1)-Y1 (J/K-1))/2.0 
30  CONTINUE 

XEr(  J/1 ) =(-3.*X1  (J,1 ) +4. *X1 (J, 2) -XI (J,3) ) /2.0 
YEF(J/1)=(-3.*Y1(J/1)t4.*Y1(J/2)-Y1(J,3))/2.0 
XEF(J/KMAX)=(3.*Xl(J/KMAX)-4.*X1(J/KM)+Xl(J/KM-1))/2.0 
YEF(J/KMAX)=(3.*Y1(J/KMAX)-4.*X1 ( J/KM)+Y1 ( J / KM-1 ) ) / 2 . 0 
40  CONTINUE 

DC  50  K=1/ KMAX/KV 
DO  50  J*1 / JM A X 

DINV=XXF(J,K)*YEF(J/K)-XE=(J/K)*YX=(J,K) 

I F ( D I N V .EQ.  0.0)  GO  TO  50 
Q ( K / 6/ J ) = 1./DINV 
50  CONTINUE 

CALCULATE  THE  PRESSURE  GRADIENTS  OVERALL  FLOW  FIELD 

IF( ITER  .GT.  1 ) GO  TO  120 
DO  30  K=1/KMAX 
DO  70  J=2/ JM 

DPDX(K/J)=0.5*(P(K/J+1)-P(K/J-1)) 

70  CONTINUE 
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DPDX(K/1)=0.5*C-3.*P(K/1)+4.*P(K/2)-P(K/3)) 

DPDX(K/JMAX)=0.5*(3.*P(K/JMAX)-4.*PCK,JM)+P(K/JM-1)) 

30  CONTINUE 

00  90  K = 1 / K M JrX 

DO  90  J=2/JM 

?XR  = P (K, J)-P (K/ J-1 ) 

PXP=P(K,J+1)-P(K,J) 

XPR=PXP*PXR 

IFCXPR  .LT.  0.0)  DPDX  CK/J)=Q. 5*(DPCX ( K/ J + 1 ) -D  P DX (K/J-1 ) ) 
C IFCXPR  .LT.  0.0)  DPDXCK, J ) =0 . 5* ( P X P-P X R ) 

90  CONTINUE 

DO  110  J=1/JMAX 
DO  100  K = 2 / K M 

DPDY(K/J)=0.5*(P(K+1/J)-PCK-1/J)) 

100  CONTINUE 

DPDY(1/J)=0.5*(-3.*P(1/J)+4.*P(2/J)-P(3/J)) 

DPDY(KMAX/J)=0.5*C3.*P(KMAX/J)-4.*P(KM/J)+P(KM-1/J>) 

110  CONTINUE 

00  1 20  J = 1 / J M A X 
DO  120  K = 2 / K M 
PYR=P(K,J)-P(K-1,J) 

PYP=P(K+1,J)-P(K,J) 

YPR=PYP*PYR 

I F ( Y ? R .LT.  0.0)  OPDYCK,J)=3.5*(OPCY(K+1/J)-DPDYCK-1/J)) 
C IFCYPR  .LT.  0.0)  DPDY(K/J)=0.5*(PYP~PYR) 

120  CONTINUE 
C 

CALAULATE  THE  TR A NS F CM AT I ON  OF  PRESSURE  GRADIENTS 
FORM  THE  WEIGHT  =UNCTION  WITH  THE  FIRST  ORDER  OF  PRESSUR 
GRADIENT 


DO  130  K=1/KMAX 
DO  130  J = 1 / J M A X 
PXD=DA3S  CDPDXU/ J)  ) 

PED  = DA3S  CDPDYCK/ J) ) 

PGX=QCK/6/J)*C3XD*YEF(J,K)-PED*YX=(J/K)) 

PQY=3CK/6rJ)*(PED*XXF(j,K)“PX0*XEF{JrK)) 

W(K,J)=DA3S(P0X)+DA3S(PQY) 

W(K/ J)=AFA  + 3TA*WU,  J) 

130  CONTINUE 

CALCULATE  THE  WEIGHTED  FUNCTION  SMOOTHLY 

DO  160  1=1/1 MAX 
00  140  K =1 /KMAX/KM 
DO  140  J = 2 / J M 

WW=0.5* ( W (X, J+1 ) +W U, J-1 ) ) 
W(K/J)=W(K/J)+0.4*(WW~W(K/J) ) 

140  CONTINUE 

W ( K / 1 )=2.*W(K/2)-WCK/3) 
W(K/JMAX)=2.*W(K,JM)-W(K/JM-1) 

DO  150  J=1 /JMAX/JM 
DO  150  K = 2 / K M 

WW=0. 5*<W(K+1 / J)+W(K-1/J) ) 
W(K/J)=W(K/J)+0.4*(WW-W(K/J)) 
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c 

c 


150 


160 


CONTINUE 


W(1/J)=2.*W(2/J)“W(3*J) 


W(KMAX/J)=2.*W(KM/J)-W(KN-1/J) 
DO  160  K=2/KM 
DO  160  J=1,JM 


WW  = 0.25*(W(K, J+1 ) +W(K,J-1 J+WCK+1, JJ+WCK-1  , 
WCK,J)=w(K,J)+0.4*(*w-W(K, 

CONTINUE 


CALCULATE  THE  DERIVATIVES  OF  WEIGHT  FUNCTION 

DO  170  K = 2^ K H 
DO  170  J = 2 ✓ J M 

WED=0.5*(W(K+1,J)-w(K-1,J)) 

WX0=0.5*(W(K, J+1 )-W(K, J-1 ) ) 

DWX(K,J)=  WXD*YEF(J,K)-WED*YXF(J/K) 

DWY(K,J)=  WED*XXF (J,X)-WXO*XEF( J,K) 

1 70  CONTINUE 
RETURN 
END 
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